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The continual development of computer technology has enabled the expansion of 
intelligent control into the field of underwater robots, where potential uses include 
oceanographic research, environmental monitoring and military mine countermeasures. 
With the naval focus shifting to operations in the littorals, and the need to lower cost of 
operations, tetherless autonomous vehicles are now being proposed for use in very 
shallow water minefield reconnaissance. These areas are dominated by a highly 
energetic environment arising from waves and currents. Motion control in such an 
environment becomes a difficult task and is the subject of this work. 

The main objective of this dissertation, is to show through modeling and 
simulation, and in-ocean experimental validation, that intervention tasks performed by 
intelligent underwater robots are improved by their ability to gather, learn and use 
information about their working environment. Using a new generalized approach to the 
modeling of underwater vehicles, which directly includes disturbance effects, a new 
Disturbance Compensation Controller (DCC) is proposed. The DCC, employing onboard 
vehicle sensors, allows the robot to learn and estimate the seaway dynamics. This 
self-derived knowledge is embedded in a non-linear sliding mode control law which 
allows significantly improved motion stabilization. The performance of the DCC has 
been verified in Monterey Harbor using the NPS Phoenix AUV. 



v 




VI 





TABLE OF CONTENTS 



I. INTRODUCTION 1 

A. MOTIVATION 1 

B. BACKGROUND 2 

C. OBJECTIVES 5 

D. CONTRIBUTIONS 5 

E. DISSERTATION ORGANIZATION 6 

II. MATHEMATICAL MODELING OF UNDERWATER VEHICLES 9 

A. INTRODUCTION 9 

B. COORDINATE SYSTEMS, POSITIONAL DEFINITIONS AND 

KINEMATICS 10 

1 . Reference Frames 11 

2. Euler Angles 13 

3. Kinematics 14 

C. NEWTON’S SECOND LAW 17 

1 . Translational Equations of Motion 18 

2. Rotational Equations of Motion 19 

3. Equations of Motion in Free Space 20 

D. RESTORING FORCES AND MOMENTS 21 

E. WAVE THEORY AND HYDRODYNAMICS 23 

1 . Linear Wave Theory 23 

2. Hydrodynamics 29 

F. HYDRODYNAMIC FORCES AND MOMENTS 32 

1 . Radiation Induced Forces and Moments 33 

2. Excitation Forces 36 

3. Viscous Damping Forces and Moments 37 

G. CONTROL FORCES 39 

1 . Propulsion Forces , 40 

2. Actuator Forces 43 

vii 



H. 6DOF EQUATIONS OF MOTION 44 

1 . No Fluid Motion 44 

2. Modifications To Account For Fluid Motion 44 

I. DEVELOPMENT OF LONGITUDINAL SURGE MODEL 47 

J. SUMMARY 48 

HI. DISTURBANCE ANALYSIS 49 

A. INTRODUCTION 49 

B. STATISTICAL DESCRIPTION OF SEA WAVES 50 

1. Wave Spectral Densities 5 1 

2. Statistical Description of Wave Amplitudes 52 

3. Empirical Formulation of Sea Spectra 54 

C. LINEAR REPRESENTATION OF SEA WAVES 58 

1 . Spectral Approximations 58 

2. State Space Formulation 59 

3. Spectral Modifications for a Moving Vehicle 60 

D. WAVE MODELING USING RECURSIVE METHODS 62 

E. APPLICATION TO DISTRIBUTED SIMULATIONS 68 

F. SUMMARY 72 

IV. PARAMETER ESTIMATION 73 

A. INTRODUCTION 73 

B. ESTIMATION METHODOLOGY 74 

C. PARAMETER IDENTIFICATION MODEL 75 

D. SYSTEM IDENTIFICATION 77 

1. Input Signal Design 77 

2. Data Collection 79 

3. Parameter Identification Results 81 

E. PROPULSION SYSTEM BALANCING 85 

F. SUMMARY 87 

V. DISTURBANCE REJECTION THEORY 89 

A. INTRODUCTION 89 

viii 



B. OVERVIEW OF CLASSICAL CONTROL TECHNIQUES 89 

1. Feedback Control 90 

2. Feedforward of Directly Measured Disturbance 91 

C. MODEL BASED CONTROL 92 

1. Linear Quadratic Regulator Control 93 

2. LQR Control with Disturbance Feedforward 94 

3 . N onlinear Methods 95 

D. CASE STUDIES IN DISTURBANCE REJECTION 98 

1. Case I. High Gain LQR Control 98 

2. Case II. LQR Control with Disturbance Estimation 

Feedforward 107 

3. Case HI. Sliding Mode Control with Measured Disturbance 

Feedforward 114 

4. Disturbance Rejection Case Comparison 1 19 

E. SUMMARY .' 122 

VI. DISTURBANCE COMPENSATION CONTROLLER (DCC) 125 

A. INTRODUCTION 125 

B. DCC OVERVIEW 125 

C. STATE AND DISTURBANCE ESTIMATION 126 

1. Model and Filter Development 127 

2. Kalman Filter Algorithm. 128 

3. Asynchronous Data Processing 129 

D. ASYNCHRONOUS SIMULATOR DEVELOPMENT 129 

E. INITIAL IN- WATER TESTING 132 

F. SUMMARY 139 

VII. ESTIMATION OF WAVE DIRECTIONALITY FROM A MOVING 

PLATFORM 141 

A. INTRODUCTION 141 

B. WAVE SPECTRA AND DIRECTIONAL ESTIMATES 14 1 

C. WAVE BUOYS 142 

ix 



D. EXTENSION TO SUBSURFACE SENSORS 145 

E. CORRECTION FOR A MOVING PLATFORM (THE DOPPLER 

EQUATION) 147 

F. DETERMINATION OF CONTROL COMMANDS 149 

G. SUMMARY 149 

VIII. EXPERIMENTAL RESULTS 15 1 

A. INTRODUCTION 151 

B. SOFTWARE IMPLEMENTATION 15 1 

C. EXPERIMENTAL VALIDATION OF THE DCC 153 

D. WAVE DIRECTIONAL ESTIMATION USING THE NPS PHOENIX 

AUV 160 

E. SUMMARY 161 

IX. CONCLUSIONS AND RECOMMENDATIONS 165 

A. CONCLUSIONS 165 

B. RECOMMENDATIONS FOR FUTURE RESEARCH 166 

APPENDIX A. EQUATIONS OF MOTION AND PARAMETERS FOR THE NPS 

PHOENIX 169 

APPENDIX B. DOPPLER SENSORS 179 

APPENDIX C. THE NPS PHOENIX AUV 187 

APPENDIX D. AUVFEST '98 201 

REFERENCES 211 

BIBLIOGRAPHY 219 

INITIAL DISTRIBUTION LIST 239 



x 



LIST OF FIGURES 



Figure 1 . 1 Phoenix AU V Control Architecture 7 

Figure 2.1 Coordinate Frames and Axes Convention 12 

Figure 2.2 Monochromatic Progressive Surface Gravity Wave [Rahman 1994] 24 

Figure 2.3 Particle Orbits And Variation Of Particle Velocity Amplitude With 

Depth [Sarpkaya 1981] 28 

Figure 2.4 Load Regimes 31 

Figure 2.5 Normalized Drag Force vs. Angle of Attack 40 

Figure 2.6 Thrust and Torque Coefficients versus Angle of Attack [Lewis 1988] 42 

Figure 3.1 Wave Spectral Density Comparisons [Berteaux 1976] 52 

Figure 3.2 Comparison of Empirical Spectra 57 

Figure 3.3 Incident Wave Directions 61 

Figure 3.4 Encounter Angle 62 

Figure 3.5 Comparison Of Innovation Covariance To AR Model Order 66 

Figure 3.6 8 th Order AR Model-Fitted To Monterey Bay Data 67 

Figure 3.7 100 th Order AR Model-Fitted To Monterey Bay Data 67 

Figure 3.8 Wave Elevation Time Series, Monterey Bay April 1998 69 

Figure 3.9 Surface Elevation PSD, O nT1 (/) 70 

Figure 3.10 Sub-surface Water Particle Velocity Series, (H=45 m, Z=22.5 m) 70 

Figure 3.11 Horizontal Water Particle Velocity PSD, O uu (/), (H=45 m, Z=22.5 m) 71 

Figure 4.1 Control Input Time Series 78 

Figure 4.2 Control Input Frequency Content 78 

Figure 4.3 Example of Measured Data, Propeller Revolutions, n , And 

Ground Velocity, u g , From The RDI 80 

Figure 4.4 Digital Voltage to Propulsion Motors 80 

Figure 4.5 Parameter Evolution 82 

Figure 4.6 Covariance Evolution 83 

Figure 4.7 Autocorrelation of Estimation Filter Residuals 84 

Figure 4.8 Comparison of Measured and Simulated Vehicle Response 85 



xi 



Figure 4.9 Nonlinear Shaft Revolution to Motor Voltage Function 86 

Figure 5.1 Block Diagram of a Feedback Controller 90 

Figure 5.2 Feedforward-Feedback Controller Block Diagram 92 

Figure 5.3 Block Diagram of State-Space System 92 

Figure 5.4 Block Diagram of Closed-Loop State-Space System with LSF 93 

Figure 5.5 Block Diagram of Closed-Loop State-Space System with LSF and 

Estimated Disturbance Feedforward 96 

Figure 5.6 Sample Disturbance Input Time Series 99 

Figure 5.7 Position Response for Case la with Monochromatic Disturbance Input 102 

Figure 5.8 Propeller Response for Case la with Monochromatic Disturbance Input .... 102 
Figure 5.9 Position Response for Case lb with PM Spectrum Based Disturbance 

Input 103 

Figure 5.10 Propeller Response for Case lb with PM Based Disturbance Input 104 

Figure 5. 1 1 Position Response for Case Ic with Monterey Bay Disturbance Input 105 

Figure 5.12 Propeller Response for Case Ic, Monterey Bay Disturbance Input 106 

Figure 5.13 Position Response for Case Ic, Monterey Bay Disturbance Input 106 

Figure 5.14 Propeller Response for Case Ic, Monterey Bay Disturbance Input, 

rpms Within Design Limits 107 

Figure 5. 15 Comparison Of Control Input Covariance To Normalized Vehicle 

Position Covariance, Monterey Bay Wave Data 108 

Figure 5.16 Position Response for Case Ha, Monochromatic Disturbance Input 109 

Figure 5.17 Propeller Response Case Ha, Monochromatic Disturbance Input 1 10 

Figure 5.18 Position Response for Case lib, PM Based Disturbance Input Ill 

Figure 5.19 Propeller Response for Case lib, PM Based Disturbance Input 112 

Figure 5.20 Position Response for Case lie, Monterey Bay Disturbance Input 1 12 

Figure 5.21 Propeller Response for Case lie, Monterey Bay Disturbance Input 113 

Figure 5.22 Position Response for Case lie, Monterey Bay Disturbance Input, 

rpms Within Design Limits 113 

Figure 5.23 Propeller Response for Case lie, Monterey Bay Disturbance Input, 

rpms Within Design Limits 1 14 

xii 



Figure 5.24 Disturbance Cancellation Case Ilia, top to bottom respectively, position 
vs. time, propeller RPM vs. time, and a phase plane plot of the 

sliding surface 116 

Figure 5.25 Controller Performance Comparison, For A Controller That Uses All 
The Disturbance Components (Dashed Line), And A Controller That 

Uses Only Fluid Velocity For Disturbance Cancellation (Solid Line) 1 17 

Figure 5.26 Position Response for Case Illb, PM Based Disturbance Input 118 

Figure 5.27 Propeller Response for Case IHb, PM Based Disturbance Input 1 19 

Figure 5.28 Position Response for Case IHc, Monterey Bay Disturbance Input 120 

Figure 5.29 Propeller Response for Case IHc, Monterey Bay Disturbance Input 120 

Figure 5.30 Comparison Of Controllers For Various Gains, PM Based Disturbance 

Input 121 

Figure 5.3 1 Vehicle Frequency Response, (Disturbance Input To Position 

Output), Superimposed Over The PM Based Disturbance Input 122 

Figure 6.1 Block Diagram of Disturbance Compensation Controller (DCC) 126 

Figure 6.2 Asynchronous simulation with realistic noise models - Disturbance from 

PM Spectrum, H s = 1 meter, operating depth 10 meters 130 

Figure 6.3 Simulated and Estimated Position Response, Using Final DCC Design 131 

Figure 6.4 Simulated and Estimated Thrust Response, Using Final DCC Design 131 

Figure 6.5 DCC Error Covariance Evolution 132 

Figure 6.6 Short Segment In-Water Results, Position for Radv=10 .’ 133 

Figure 6.7 Short Segment In-Water Results, Relative Velocity for Radv=10 133 

Figure 6.8 Short Segment In-Water Results, Fluid Velocity Estimate for Radv=10 134 

Figure 6.9 Short Segment In-Water Results, Thrust Estimate for Radv=10 134 

Figure 6.10 Short Segment In-Water Results, Propeller RPMs for Radv=10 135 

Figure 6. 1 1 DCC Frequency Response, ADV Input to Propeller Output , Radv= 10 .... 136 
Figure 6.12 DCC Frequency Response, ADV Input to Propeller Output , Radv=100.. 136 

Figure 6.13 Short Segment In-Water Results, Position for Radv=100 137 

Figure 6.14 Short Segment In-Water Results, Relative Velocity for Radv=100 137 

xiii 



Figure 6.15 Short Segment In-Water Results, Fluid Velocity Estimate for 

Radv = 100 138 

Figure 6.16 Short Segment In-Water Results, Thrust Estimate for Radv=100 138 

Figure 6.17 Short Segment In-Water Results, Propeller RPMs for Radv=100 139 

Figure 8.1 Software Implementation of DCC 152 

Figure 8.2 Block Diagram of DCC Implementation 152 

Figure 8.3 Comparison of Measured and Estimated Position, April 22, 1999, 

Run # 3 154 

Figure 8.4 Propeller Response, April 22, 1999, Run # 3 155 

Figure 8.5 Measured Ground Velocity, April 22, 1999, Run # 3 155 

Figure 8.6 Comparison of Measured and Estimated Relative Velocity, 

April 22, 1999, Run # 3 156 

Figure 8.7 Fluid Velocity Estimate, April 22, 1999, Run # 3 156 

Figure 8.8 Estimated Thrust, April 22, 1999, Run # 3 157 

Figure 8.9 Comparison of DCC Performance, Simulation and Experimental 158 

Figure 8.10 Comparison of Transient Response for Various Control Gains 159 

Figure 8. 1 1 Transient Response Prediction of the DCC 160 

Figure 8.12 Sample Direction Energy Plot From Phoenix Data, Nov. 4, 1998 

(Run # 9), Gulf of Mexico 162 

Figure 8.13 Sample Direction Spectrum Plots From Phoenix Data, 

Nov. 4, 1998 (Run # 9), Gulf of Mexico 162 

Figure 8.14 Short-term Current Estimation from Phoenix, November 8, 1999, 

(Run # 2), Gulf of Mexico 163 

Figure B.l ADVOcean Probe 179 

Figure B.2 ADV Sampling Volume 180 

Figure B.3 Standard ADVOcean Probe 181 

Figure B.4 ADVOcean Probe with Optional Sensors Housing 182 

Figure B.5 ADVOcean Dimensions 184 

Figure C.l Physical Layout of Phoenix AUV 187 

Figure C.la Wire Frame Diagram of New NPS AUV 188 



xiv 



Figure C.lb Solid Model of New NPS AUV 188 

Figure C.2 Brushless DC Motors 191 

Figure C.3 Old 3 -in. Props 191 

Figure C.4 Present 7-in Ducted Props 192 

Figure C.4a Present 3.5-in Thrusters 192 

Figure C.5 Phoenix Control Actuators 193 

Figure C.6 Mission Control Computer 193 

Figure C.7 FreeWave Modem : 194 

Figure C.8 Systron Donner MotionPak 194 

Figure C.9 PrecisionNav TCM2 Compass 195 

Figure C. 10 RDI Navigator DVL 195 

Figure C.ll ADVOcean Acoustic Doppler Velocimeter 196 

Figure C.12 ST-725 and ST-1000 Sonars 196 

Figure C.13 Transport Cart 197 

Figure C.14 Mobile Lab Internals 198 

Figure C.15 Mobile Lab 199 

Figure C.16 Support Craft 199 

Figure C.17 LXT Acoustic Tracking System 200 

Figure D.l Phoenix Shipping Containers 202 

Figure D.2 Packaging of Support Equipment 203 

Figure D.3 Gulfport Temporary Work Space 203 

Figure D.4 R/V Gyre at Anchor 204 

Figure D.5 Phoenix Being Launched from the Gyre 204 

Figure D.6 AUVFEST Operations Area 205 

Figure D.7 Clusters Identified 209 

Figure D.8 Lower right cluster with centroid 210 



xv 





XVI 






LIST OF TABLES 



Table 2. 1 Standard Underwater Vehicle Notation 1 1 

Table 2.2 Summary Of Small Amplitude Linear Wave Equations 27 

Table 3. 1 Wave Amplitude Means 53 

Table 3.2 Expected Maximum Amplitudes [Longuet-Higgins 1953] 54 

Table 3.3 Sea State and Wave Parameter Comparison [Berteaux 1976] 55 

Table 4. 1 Parameter Estimation Values and Residual Statistics 82 

Table 8. 1 Sample Summary of DCC Validation Runs 153 

Table B.l Standard ADVOcean Features 182 

Table B.2 ADVOcean Performance Specifications 183 

Table B.3 Navigator DVL Specifications 185 

Table C.l Vehicle Upgrades and Acquisitions 189 

Table C.l Vehicle Upgrades and Acquisitions (Cont.) 190 

Table D. 1 Sample Phoenix Missions 206 



xvii 



XV111 



ACKNOWLEDGEMENTS 



I would like to express my most sincere gratitude to all the members of my 
Doctoral Committee for their constructive guidance and suggestions during the course of 
this research. In particular, I wish to thank Professor Anthony Healey, my dissertation 
supervisor, for his time, patience, and enthusiasm, without which, this work would not 
have been possible. His vast knowledge, experience, and professionalism have been an 
inspiration. 

A special thanks to Dr. David Marco for his advice, guidance, and programming 
pertaining to the real-time implementation of the control code in this research. Without 
your time and effort I would never have made it. 

I would also like to recognize the financial support of the Office of Naval 
Research (Mr. Tom Curtin) under contract No. N0001498WR30175. 

Finally, I dedicate this dissertation with love and thanks to my wife Ellen and our 
two beautiful children Taylor and Chelsea. Your understanding and support during the 
frequent trips and long hours have made you a part of this work. 



xix 



I. INTRODUCTION 



A. MOTIVATION 

The continual development of computer technology has enabled the expansion of 
intelligent control into the field of remotely operated underwater vehicles. With this 
advent, tetherless Autonomous Underwater Vehicles (AUVs), have arrived and are 
expected to be self-sufficient with respect to power and control. Driven by the need to 
lower the cost of operations, applications associated with these vehicles are both 
appealing and numerous. The potential uses for these vehicles include but are not limited 
to: scientific (oceanography, geology, geophysics,... [Curtin 1993, Smith 1994 and 
Pereira 1996]), environmental (waste disposal monitoring, wetland 
surveillance,... [Chase 1998]), commercial (oil and gas, submerged cables, 
harbors,... [Hartley 1991, Butler 1998, Kato 1998]), or military (mine warfare, tactical 
information gathering, smart weapons,... [Honegger 1996]). During a mission, an AUV 
is expected to carry sensors (side scan sonar, bathymetry, bottom profiler,), track to a 
certain planned trajectory (traveling from point A to point B, performing the systematic 
exploration of a zone,...), and even make on-line decisions allowing for mission 
reconfiguration, [Yavani 1996, Byrnes 1996]. 

The payload, which the vehicle carries, can place constraints on the motion of the 
vehicle and mission planning. With new emphasis on naval mine countermeasures and 
reconnaissance in very shallow waters (VSW) [ONR 1998], a typical mission may 
require a vehicle to follow an unknown and pro file- varying bottom at a desired altitude, 
while maintaining a nominal constant forward speed, and avoiding obstacles when 
necessary, all while operating in a highly energetic environment. Since small AUVs are 
sensitive to wave induced motions [An 1996], AUVs in VSW will perform better if they 
derive knowledge, using their onboard sensors, about their environment and operating 
area. 

The question then is: How can environmental knowledge be obtained and used in 
a real-time control architecture in order to correctly and safely achieve an assigned 
mission? The integration of sensory data in closed-loop control of mobile robots or 
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vehicles has been a topic of recent research in both control architecture [Saridis 1983, 
Healey 1996, Marco 1996] and execution [Marco 1998] . Based on this work, two broad 
categories can be distinguished: 1) the behavioral approach and 2) the sensor based 
control approach. The first is too restrictive for general use, since it needs to specify the 
problem as a set of simple input/output relations for all possible situations. Among the 
second class, there are several methods proposed, including intelligent servo control 
[Marco 1996c], path following [Fryxell 1995] and collision avoidance [Williams 1990], 
and the task function approach [Santos 1995]. 

The task function approach allows a framework by which various control laws 
may be develop for specific mission tasks. The embedded control laws are then used 
based on sensor input decisions with a mission manager controlling the commands and 
task assignments. The goals of this research include the need to show that task 
assignments can be performed with less environmentally induced motions if the 
characteristics of the operating environment are known. This is particularly important 
because of the need to operate in environmentally energetic areas and because of the 
unpredictable nature of the VSW region. 

Classically, this problem falls into the broad category of servo control disturbance 
rejection. This problem is well known and has a rich history of study. In principle, a 
high gain feedback servo is known to track commands and reject disturbances (unwanted 
inputs). Requirements for loop stability and sensor noise rejection are competing and the 
subject of optimal control formulations. It is conjectured that knowledge or direct 
measurement/estimation of these disturbances improves control performance. 
Measurement/estimation of ocean disturbances is both difficult, as well as subject to the 
unpredictable nature of the ocean. For all the ongoing work in the offshore oil industry, 
wave disturbances are still not measured and used directly for control purposes. 

B. BACKGROUND 

Autonomous systems for small unmanned untethered underwater vehicles have 
several features that separate them from traditional marine control and guidance systems. 
The single most important feature is that it is desirable to control a small underwater 
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vehicle with relatively high velocities along or about two or more axes at the same time 
reliably. This leads to stronger coupling, larger nonlinearities, more state equations and 
more unknown parameters in the vehicle's equations of motion than what would be 
present with surface vessels. This single most important fact makes the control of small 
underwater vehicles relatively difficult. 

Other factors include but are not limited to: 

• A small underwater vehicle may be controllable in all six degrees of freedom, 

• A small underwater vehicle has a natural frequency in the same range as the 
environmental disturbances; 

• Actuator dynamics are much faster on small underwater vehicles; 

• Power for control and operations is limited by what the vehicle can carry 
onboard; 

• Man- in-the- loop operations and human intervention for fault response are not 
possible, if controller problems develop; 

• Small underwater vehicles, having a higher bandwidth, may compensate for 
the first order wave disturbances by means of feedforward and feedback 
control laws. 

There are numerous factors that must be considered in the design of a control 
system that is capable of being implemented in an operating vehicle. Some of the most 
important items are: 

• The bandwidth is limited by the control system’s sampling rate, computational 
time delays and sensor communication delays. The time delay factor is 
extremely important when the control law is dependent on information from a 
hydro-acoustic source. 

• The control system's sensitivity to noise has to be considered, because it 
determines how good a sensor must be and how robust the control algorithm 
is with respect to noisy measurements. 

• The required accuracy of the system is usually a function of the commanded 
tasks; hence it is not meaningful to discuss this factor before a task 
specification is made. 
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• Stability is of the utmost importance. Depending on the type of system (SISO 
or MIMO) and the type of controller to be employed (linear or nonlinear), 
various techniques are available to rate these criteria. 

• Requirements to handle system failures, i.e., sensor failure, actuator failure 
and even computer failure, must be addressed. 

• Parameter variation is also an important aspect in robustness analysis of a 
control system. Parameter variations can be handled by a well-designed 
adaptive algorithm or by employing a control scheme that is robust to these 
variations. 

The first modem DP system was first introduced by [Balchen et al., 1976], where 
a Kalman filter approach for solving the wave noise-filtering problem was employed. 
Before the employment of this system, DP systems were mainly based on PID-controllers 
and matching notch filters which filtered out the wave "noise" on the sensors. Kalman 
filter based DP systems now dominate the market and the literature [Balchen 1976, 
Jenssen 1980, Saelid 1983, Fung 1983, Triantafyllou 1983, Fossen 1995, and Sorensen 
1996]. 

The Kalman filter based DP systems consist of one high and one low-frequency 
model where the control action is based on the low frequency model. The high frequency 
model is used to filter out the relatively high-frequency, first-order, wave noise that 
appears on the sensors. 

This filtering is not necessary for small underwater vehicles, since in theory, 
compensation for the first order wave disturbance’s is possible. For small vehicles, the 
DP case becomes a special case of the general control problem where desired trajectory 
velocities are zero. Although there exist many available DP systems for surface vessels, 
the Simrad-Albatross as one example, and even several ROV DP systems, such as the 
[Marquest 1991] system and the systems developed and deployed on the Norwegian 
Experimental Remotely Operated Vehicle (NEROV) by [Fossen and Satagun 1991], 
there currently are NO KNOWN AUTOMATIC DP systems for AUVs. 
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C. OBJECTIVES 



The main objective of this research is to show through modeling, simulation and 
experimental validation that intervention tasks performed by intelligent underwater 
robots are improved by their ability to gather, learn and use information about their 
working environment. 

Specifically this dissertation addresses the following topics: 

• Generation and verification of mathematical models; 

• Measurement and use of wave disturbance information for control 
compensation; and 

• Implementation and verification of real-time embedded control processes. 

D. CONTRIBUTIONS 

This dissertation contains both theoretical and experimental contributions to the 
field of applied underwater vehicle control. The theoretical contributions are partly 
found in the modeling chapter (Chapter II), the disturbance rejection chapter (Chapter V), 
and the wave direction estimation chapter (Chapter VII). The experimental contributions 
(Chapter VIII) are from work carried out during operational testing of the NPS Phoenix 
AUV. The Phoenix and its subsystems, including system upgrades necessary to perform 
this research, are described in Appendix C and in [Marco 1996]. 

The contributions in this dissertation are described below: 

• Chapter II provides a new generalized approach to the modeling of small 
underwater vehicles subject to shallow water wave and current effects. Using 
appropriate modeling formulations, as opposed to adding white or colored 
noise, random disturbances, or "RAO" based motions, the fluid effects are 
incorporated directly into the equations of motion. In addition, this 
formulation provides the ability to construct a generalized distributed 
simulation capability for the evaluation of underwater vehicle systems in 
shallow water. (Generally useful to U.S. Navy tactical decision making). 

• It is proven through simulation (Chapter V), and verified by experimental 
validation (Chapter VIII), that measuring the water column velocities directly, 
wave induced disturbances can be substantially compensated by the vehicle 
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control system. This technique eliminates the need to develop and incorporate 
sophisticated predictive disturbance models in the control system design. 

• In Chapter VII, it is shown how small underwater vehicles, using direct fluid 
measurements, can obtain short-term wave magnitude, directionality and 
current estimates, thereby providing useful information in the area of tactical 
oceanography. It is also shown how a general seaway direction may be 
obtained for use in mission planning and control. 

E. DISSERTATION ORGANIZATION 

This dissertation is aimed at contributions in the area of control of intelligent 
underwater robots. Figure 1.1 represents the general control architecture for the Phoenix 
AUV. Using this figure as a roadmap, the chapters and appendices consider the 
following: 

• Chapter II. "Mathematical Modeling Of Underwater Vehicles" derives the 
complete set of nonlinear equations of motion for a small underwater vehicle 
subject to shallow water waves. Kinematics, Newton's laws of angular and 
linear momentum, general hydrodynamics and external force modeling are 
discussed in detail. 

• Chapter III. "Disturbance Analysis" describes how deterministic and 
stochastic disturbances can be modeled for use in the vehicle control 
architecture. Statistical description, state space representation and 
autoregressive (AR) modeling are used to illustrate these ideas. 

• Chapter IV. "Parameter Identification" describes the theoretical foundation 
and experimental results associated with determining the parameters and 
coefficients used to model the NPS Phoenix AUV in the surge direction. 
Comparisons between vehicle experiments and simulations show the level of 
certainty associated with the identified parameters. 

• Chapter V. "Disturbance Rejection Theory" provides a survey of the 
classical and modem approaches to disturbance rejection and compensation. 
This chapter describes how a nonlinear surge controller can be applied to an 
underwater vehicle design. It details three specific controller designs: a LQR 
design, a LQR with embedded disturbance model design, and a nonlinear 
sliding mode controller with disturbance compensation. 
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Figure 1.1 Phoenix AUV Control Architecture 



• Chapter VI. "Disturbance Compensation Controller (DCC)" is a discussion 
of the design and implementation of a combined nonlinear model-based 
extended Kalman filter with a SMC for dynamic positioning. The state 
estimator is critical to the controller performance since the surge controller 
requires states that are not measured with existing systems and sensors. This 
chapter compares actual signals with estimated signals, and shows that by 
properly designing the filter gains, the measurement noise levels are reduced 
and the controller bandwidth is increased. Also, critical to the point of the 
dissertation, ocean fluid particle velocities, the disturbance, are estimated. 

• Chapter VII. "Wave Directional Estimation From A Moving Platform" 
outlines the theory used to estimate wave directions from an AUV. It also 
describes the method by which a heading command, based on the dominant 
wave energy direction, is calculated. This heading command is used in the 
heading controller to properly orient the vehicle in the direction of the 
maximum energy, thereby reducing the drag force on the vehicle. 

• Chapter VIII. "Experimental Results and Validation" is an analysis of the 
experimental missions performed during this research. It describes the 
real-time implementation of the various processes (filters, controllers and data 
acquisition code) operating in the Phoenix. The chapter shows that the 
Disturbance Compensation Controller (DCC) is capable of dynamic station 
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keeping in the presence of waves. It also displays and analyzes directional 
wave spectrum estimates obtained during AUVFEST '98. 

• Chapter IX. "Conclusions and Recommendations" contains comments, 
conclusions and recommendations for future work. 

• Appendix A. "Equations of Motion and Parameters for the NPS Phoenix 
AUV" provides the physical characteristics, the hydrodynamic coefficients 
and the equations of motion for the NPS Phoenix AUV. 

• Appendix B. "Doppler Sensors" outlines the specifications and principles of 
operation of the SonTek® ADV and the RDI® WORKHORSE DVL. 

• Appendix C. "The Phoenix AUV" describes in detail the vehicle more 
closely, including upgrades and acquisitions that allowed the Phoenix to 
transition from a test tank vehicle to open ocean operations. 

• Appendix D. "AUVFEST ’98" describes the NAVO sponsored underwater 
vehicle demonstration exercise that took place during November 1998 off the 
coast of Gulfport, MS. 
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II. MATHEMATICAL MODELING OF UNDERWATER VEHICLES 



A. INTRODUCTION 

The mathematical modeling of underwater vehicles can be found throughout the 
literature. The standard equations of motion for submarine simulation from the David 
Taylor Model Basin Naval Ship Research and Development Center (DTMSRDC) are 
described in detail by [Gertler 1967] and [Feldman 1979]. In [Kalske 1989] a survey of 
the motion dynamics of underwater vehicles, including ROVs and submarines, is given. 
A description of the linear and nonlinear motion dynamics of the University of New 
Hampshire Experimental Underwater Vehicle (UNH-EAVE) is found in 
[Humphreys 1982]. 

Many other papers and books have been written on this subject (see Yuh 1990, 
Healey 1992 and 1993, and Fossen 1994 for further examples). Each model has 
developed in complexity and accuracy, but each model has vehicle specific 
simplifications making it necessary to revisit this topic. Presently, there is no single 
model that combines all aspects, including environmental disturbances and their effects, 
into one generalized model for shallow water operations of small underwater vehicles. 

Accurate modeling allows for the development of precision control. Precise 
control is needed for many maneuvers and tasks, such as autonomous docking and 
recovery, target detection, identification and localization, as well as station-keeping. 

In this chapter the generalized six-degree of freedom (6DOF) equations of motion 
(EOM) for a small underwater vehicle operating in shallow water will be developed. As 
with all previous modeling work in this area, the underlying assumptions are that: 

1) The vehicle behaves as a rigid body; 

2) The earth's rotation is negligible as far as acceleration components of the 
center of mass are concerned; 

3) The primary forces that act on the vehicle have inertial, gravitational, 
hydrostatic and hydrodynamic origins, and 

4) The hydrodynamic coefficients or parameters are constant. 
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The chapter will begin by outlining the coordinate frames and the kinematic and 
dynamic relationships used in modeling a vehicle operating in free space. Next a 
discussion of linear wave theory and basic hydrodynamics will be presented. This 
discussion will set the foundation for the various force and moment expressions 
representing the vehicle’s interaction with its fluid environment. The control forces, 
resulting from propellers and thrusters and from control surfaces or fms, that enable the 
vehicle to maneuver, will be then be detailed. With the hydrodynamic and control force 
and moment analysis complete, the full six degree of freedom equations of motion will be 
formed taking into account the necessary modifications for current and wave effects. 
While reduced order models representing flight control modes are detailed elsewhere, the 
chapter will conclude with a development and discussion of the one degree of freedom 
(1DOF) surge model. This model will be used in later chapters as the basis for a 
controller that will allow an AUV to station-keep in the presence of waves. 

B. COORDINATE SYSTEMS, POSITIONAL DEFINITIONS AND 
KINEMATICS 

For underwater vehicles that operate in three dimensional space and time, it is 
necessary to describe position/orientation and velocity/rotation rate of the vehicle by six 
independent coordinates or degrees of freedom. Typically these coordinates are chosen to 
correspond to the position and orientation and their time derivatives with respect to some 
set of mutually orthogonal coordinate axes fixed to an arbitrary origin which defines a 
reference frame. Likewise, the forces/moments acting on or produced by the vehicle can 
be referenced to a set of coordinate axes. In this dissertation, standard notation, 
[SNAME 1950], is used to describe the aforementioned 6 DOF quantities and is 
summarized in Table 2.1. 

Note that by convention for underwater vehicles, the positive x-direction is taken 
as forward, the positive y-direction is taken to the right, the positive z-direction is taken 
as down, and the right hand rule applies for angles. It is convenient to group the linear 
and angular body fixed velocities into a vector quantity x, where x = [ u , v, w, p , q, r] T , 
and the global positions and Euler angles as a vector quantity z, where 
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z = [X, Y, Z, <p, 0, y/\ l . The externally applied forces and moments are grouped into a 
vector quantity /, where/ = [Xf, Yf, Zf, Kf, Mf, Nf] T . 



DOF 


Motions 


Forces and 
Moments 


Linear and 
Angular 
Velocities 


Positions and 
Euler Angles 


1 


surge 


*/ 


u 


X 


2 


sway 


Xf 


V 


Y 


3 


heave 


Zf 


w 


Z 


4 


roll 


Kf 


p 


0 


5 


pitch 


M f 


q 


6 


6 


yaw 


N f 


r 


¥ 



Table 2.1 Standard Underwater Vehicle Notation 



1. Reference Frames 

As discussed earlier and outlined in Table 2.1, to properly describe or model the 
motion of a rigid body three independent positions and angles are required, and for 
convenience, three orthogonal coordinate frames are used. First, a global frame OXYZ, as 
shown in Figure 2.1, is defined with origin O, and a set of axes aligned with directions 
North, East and Down. This produces a right-hand reference frame with unit vectors I, J, 
and K. Ignoring the earth's rotation rate in comparison to the angular rates produced by 
the vehicle's motion, it can be said that the OXYZ coordinate frame is an inertial reference 
frame in which Newton's Laws of Motion will be valid. A vehicle's position in this 
global frame will have the vector components, ro = [XI + YJ + ZK] 

Secondly, a body fixed frame of reference Oxyz, with the origin O, and unit 
vectors i, j, k, located on the vehicle centerline, moving and rotating with the vehicle is 
defined. The origin O', will be the point about which all vehicle body force will be 
computed in later sections of this chapter. The vehicle's center of gravity (mass), CG, 
and center of buoyancy, CB, do not generally lie at the origin of the body fixed frame, nor 
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are they collocated. The positional vectors of the CG and CB relative to the origin of the 
body fixed frame are Pg and pn, respectively, and can be represented in component form 
as [ xgi + ycj + ZGk] and [xsi + yBj + ZbK\. The location of the center of mass (gravity) is 
important because Newton's Laws of Motion simplify when forces and moments on a 
body are applied to its center of mass. The center of buoyancy is determined by the 
shape of the submerged portion of the body, while the center of gravity is determined by 
the distribution of the weight. 

Lastly, a fluid frame is defined with origin Of. This frame is aligned parallel to 
the global reference frame but moves with the local fluid particles. Defining the fluid 
reference frame in this manner simplifies the hydrodynamic force modeling which will be 
discussed in later sections. 




Figure 2. 1 Coordinate Frames and Axes Convention 
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2. Euler Angles 

The transformation from one Cartesian coordinate system to another can be 
performed by means of three successive rotations in sequence. While there are several 
different methods to describe the attitude of a vehicle in the global reference frame, the 
most common method uses so-called Euler angles. This method uniquely defines the 
angular orientation of the vehicle reference frame relative to the global reference frame. 
There are several Euler angle conventions to describe these three successive rotations, see 
[Craig 1989] for a complete list. However, for the purpose of this dissertation, the so 
called "roll, pitch, yaw,” "X-Y-Z fixed angle" or "Tait-Bryant" convention will be used. 
One concern that may arise with the use of Euler angles is that a singularity exists when 
the elevation reaches 90 degrees. This limitation which can sometimes, although rarely, 
cause trouble in flight simulations and control computations may be overcome by the use 
of quaternions [Craig 1989], which introduce four rather than three variables to describe 
angular position. 

For the "roll, pitch, yaw" convention, a forward transformation is performed by 
beginning with a vector quantity originally referenced in the global reference frame. 
Then, through a sequence of three rotations it is transformed into a frame that is assumed 
to be parallel to and moving with the vehicle body coordinate frame. To start the 
transformation, begin by defining an azimuth rotation y/, as a positive rotation about the 
global Z-axis. Next define a subsequent rotation 6, (positive up) about the new 7-axis, 
followed by a positive rotation 0, about the new X-axis. The triple rotational 
transformation in terms of these three angles, is then sufficient to describe the angular 
orientation of the vehicle at any time. 

It follows that any position vector, r Q , in an original global reference frame given 
by r Q = [Xo, Y 0 ,Z 0 ], will have different coordinates in a rotated frame when an azimuth 
rotation by angle y/, is made about the global Z-axis. 

If the new position is defined by r; = [X lt Yj ,Zj], it can be seen that there is a 
relation between the vector's coordinates in the new reference frame with those that it had 
in the old reference frame. It follows that, 
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X, = X 0 cos^ + F 0 sin^ (2.1) 

F, = -X 0 sin y/ + Y g cos y / , (2.2) 

with Z/=Z 0 . This relation can be expressed in matrix form by the rotation matrix 
operation, 

rj = [T Vt z ] r 0 ; (2.3) 



where the rotation matrix T ¥ z, represents an orthogonal transformation. Multiplication 
of this rotation matrix with any vector, r 0 , will produce the components of the same vector 
in the rotated coordinate frame. Continuing with the series of rotations results in a 



combined total rotational transformation. 



T 1 (<t>,d,y) = T(<P)T(d)T(y). 



In expanded notation Equation 2.4 takes the form 



T x {(p,6,y/) = 



cy/cd 


sy/cd 


-s$ 


cy/sOstp- sy/ctp 


sy/sds(p + cy/c(p 


cOstp 


clf/sdc(p + sy/s<p 


s\f/sdc<p-cy/s<f) 


cdap 



(2.4) 



where c and 5 are abbreviations for cos and sin. It can be said that any positional vector 
in an original reference frame may be expressed in a rotated frame with coordinates given 
by the operation, 

fijk = Ti((f), 6, if/) rjjK. (2.5) 

3. Kinematics 



Kinematics deal with the relationships of motion quantities regardless of the 
forces induced by their prescribed motions. Description of a body's position both 
translational and rotational, will need to be related to velocities, both translational and 
rotational, prior to extending the analysis to accelerations. The connection between 
translational velocity and the rate of change of translational position is straight forward. 

Define a global velocity as, 



r = 



X 

Y 

Z 



( 2 . 6 ) 
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This vector will have components that are different when seen in a body fixed frame. 
Selecting the linear components of the body fixed velocity vector, v = [u, v, w] T , these 
components, in global quantities are found, using the forward transformation defined in 
Equation 2.5, to be 



u 




~X 


V 


— T x (<p,6,y/) 


Y 


w 




Z 



(2.7) 



It is a simple reverse coordinate transformation from body fixed to global coordinates to 
see that. 



~X 




u 


Y 


~T\~' (<p,d,y/) 


V 


Z 




w 



( 2 . 8 ) 



where J, 1 (0, 6, y/) is the transpose of Equation 2.4, since T x {(j),d,y /) is orthogonal. This 



shows that the progression of a vehicle in the global frame clearly depends on its local 
velocity components and its attitude. 

The connection between angular attitude and angular velocity is not quite so 
obvious. At first sight, it is tempting to define the instantaneous angular velocity of the 
vehicle simply as the rate of change of its angular position defined by the Euler angles. 
This is erroneous however, because the rotation 6, was defined as a rotation about the 
intermediate frame after a rotation y/, had been made. It should be noted that the order of 
multiplication in Equation 2.4 must be followed since the rotations do not commute, i.e., 
T(<p,6,y/)*T(6,y/,(f>) . Since the rotation do not commute, the set cannot form a vector 
space and therefore the derivatives cannot represent the body's angular velocity. 

Vehicle inertial angular rates are defined in terms of components that have 
angular velocities about the global axes. It is necessary to relate both Euler angles and 
their rates of change, to angular velocity components about the global axes and to their 
components lying along the body fixed axes in any attitude. The prime reason for this is 
that it is difficult to construct physical sensors to measure rates of change of Euler angles. 
(Rate gyros in common use today are easily constructed to measure the components of 
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the inertial angular velocity of a vehicle that lie along the vehicle's body axes.) It follows 
that the instantaneous angular velocity of the vehicle can be related to the instantaneous 
rate of change of angular orientation only after considerations of the intermediate 
transformations used. In other words, if a is defined as the angular attitude vector, 
a-[<p,0, y/f , and the inertial components of the vehicle angular rate lying along the body 
axes are defined as co= [p,q,r] T , then a = 



The details of the nonlinear functional relations involved are provided by viewing 
the rate of change of the rotation y/, as a vector quantity lying along the original Z-axis. 
The rate of change of the angle #, is viewed as a vector quantity lying along the 7-axis of 
the first intermediate frame, and the rate of change of the angle <p, is viewed as a vector 
lying along the Z-axis of the final (body fixed) frame. Each of these component rates of 
change of angular position has component parts that project onto the final frame and it is 
the sum total of all the components that give the total angular velocity as seen in the final 
frame of reference. Using the required transformations for the rate components from 
each Euler angle, 
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(2.9) 



the body rates are obtained with the result. 
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cos# 


sin# 
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Inverting Equation 2. 10 yields a solution for the rates of change of the Euler angles in 
terms of the body fixed components of the angular velocity vector, 



y 




/> + gsin0tan#-i-rcos0 tan# 
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sin <p tan# cos^tan# 
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qcos<f> -rsm<p 
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cos (p - sin <p 
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(gsin 0 + rc os 0)/ cos# 
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-sin (p! cos# cos^/cos# 


r 



( 2 . 11 ) 

Notice that for small angular rotations, as expected, 
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</> = p\ 

6 = q\ 

V = r, 

as expected. As a note, it should be pointed out that unlike the transformation matrix T j, 
T 2 is not orthogonal, therefore, T 2' 1 * T? ■ 

At this point the kinematic relationships between velocities, as seen in the body 
fixed frame, and the rates of change of global positions and Euler angles have been 
defined. The resulting set of differential equations forms a consistent set in that given a 
set of vehicle velocity data versus time, its position and attitude may be computed. 
Expending and combining equations 2.8 and 2.11, this set of differential equations is 
summarized here; 



~x 




u cos (? sin ^ + v(- cos 0 sin ^ + sin 0 sin 0 cos + w(sin 0 sin ^ + cos <f> sin 6 cos y/ 


Y 




u cos 6 sin y/ + v(cos (f> cos y/ + sin 0 sin 0 sin y/) + w(- sin (f> cos if/ + cos (p sin 9 sin y/) 


Z 




— u sin<? + vsin0cos0 + wcos0cos0 


4> 




p + q sin 0 tan 6 + r cos <p tan 6 


6 




qc os^-rsin0 


k 




(qsuup + rcos<p)/cos0 



In vector-matrix form, this set of equations can be expressed as; 




C. NEWTON’S SECOND LAW 

Since the motion of the vehicle is composed of both translational and rotational 
components, it is necessary to retain the distinction between the vehicle or body fixed 
frame and the global or inertial frame. This is particularly important because the 
dominant forces that act on a submerged body depend on the local motion of the vehicle 
relative to the fluid particles, not on the global velocities. Returning to Figure 2.1, the 
global position, velocity and acceleration vectors of the vehicle’s center of gravity is 

dr r d 2 r r 

defmed as r G = r Q .+p G , — — and — respectively. Since the center of mass lies in a 

dt dt l 
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body that is both translating and rotating, the total time derivative of r G comes from two 
parts. The first part is from the time derivative as if the body fixed frame was not 
rotating, while the second part is due to the rotation of the body fixed frame. A detailed 
explanation of this can be found in [Greenwood 1988]. 

For a general vector r, in a frame rotating at an angular velocity ca, the total 
derivative is given as 

dr 

— = r + axr, ( 2 . 12 ) 

dt 

where the first right hand side term is a rate of growth part and the second right hand side 
term is a rate of transport part. Using this logic, the velocity of the vehicle’s center of 
mass is given by 

dr r 

-£- = r 0 .+coxp c =v 0 , + coxp c . (2.13) 

dt 



The inertial velocity of the origin of the body fixed frame is expressed as vo- and may be 
written in either global or body fixed coordinates, therefore, 

dX dY dZ 

v 0 .=r 0 .=[—I+—J+—K] = [ui+vj+wk ]. (2.14) 

dt dt dt 

1. Translational Equations of Motion 



With the coordinate frames defined, Newton's second law of motion, 

f Trans = ~( mv c)’ ^ be used to formulate a translational model of the system. The 
dt 

global acceleration of the center of mass is derived by differentiating the velocity vector, 
r G , realizing that the center of mass lies in a rotating reference frame. Considering the 
total differential, the global acceleration of the center of mass becomes, 

r G =v 0 . + d)xp G +Q)XO)xp G +Q)Xv 0 .. (2.15) 

The translational equation of motion is found by equating the product of this 
acceleration and the vehicle mass, to the net sum of all forces acting on the vehicle in 
three translational degrees of freedom (X, Y, Z). One important factor to recognize is that 
the equation of motion derived in this manner is a vector equation with the components 
expressed in the body fixed frame with unit vectors i, j and k . As discussed earlier, this 
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has been deliberately done because the dominant forces acting on a submerged body in 
motion are developed in relation to the shape of the vehicle and are more conveniently 
expressed in relation to the body axes. Applying Newton's second law results in the 
vector equation, 

f Trans =m{x 0 ,+a)xp G +coxcoxp G +coxx 0 ,}. (2.16) 

2. Rotational Equations of Motion > 

To develop the rotational equations of motion, the sum of applied moments about 
the vehicle's center of mass is equated to the rate of change of angular momentum of the 
vehicle about its center of mass. In the practical case of marine vehicles, however, the 
statement just made is modified slightly because it is much more difficult to assess the 
vehicle's mass moments of inertia about its center of gravity ( CG ), as the CG changes 
with loading. It becomes simpler to evaluate the mass moments of inertia about the body 
fixed frame which tends to lie along the vehicle's axes of symmetry. The inertia tensor 
required to be computed is, 

h = 

where, 

A 

/„ = ^ =l dm i (y 2 +z 2 ) and I ^ = I yx = dm i (xy ) , for example. 

1=1 

The angular momentum of the body is given by, 

h 0 =I 0 6). (2.18) 

The total applied rotational moments about the vehicle's reference frame origin is given 
by ’ 

m & =h 0 .+p G x{mv G ). (2.19) 

Differentiating Equation 2.18, the rate of change of angular momentum is found to be, 

h o ,=I o . 6 )+( 0 xh o ., (2.20) 

and the acceleration of the global position vector /y, is given by, 



^ xx ^ xy 1 xz 
I yx I yy ^ yz 
I zx I zy 1 zz 



(2.17) 
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r, =v o ,+Q)Xv 0 .. (2.21) 

Substituting equations 2.20 and 2.21 into Equation 2.19, the rotational equation of motion 
in vector form is given by, 

m r 0 , =I o 0)+0)X(I o 0)) + m{p G xv o , + p G X0)Xv o .}. (2.22) 

(Again, for a more detailed discussion, the reader is referred to [Greenwood 1988]). 

3. Equations of Motion in Free Space 

At this point, there are three translational equations obtained from Equation 2.16, 
three rotational equations obtained from Equation 2.22, and six unknown velocities 
(v and c 0 ) . This set of equation written in long form is; 

m[u -vr + wq — x G ( q 2 +r 2 ) + y G ( pq -r) + z c ( pr + q )] = X f , 

m[v + ur - wp + x G (pq + r) - y G (p 2 + r 2 ) + z G (qr - p)] = Y f , 

m[w -uq + vp + x G (pr -q) + y G (qr + p)-z G (p 2 +q 2 )] = Z f , 

I x p + (I Z -I y )qr + I xy (pr-q)-I yz (q 2 - r 2 )- I x£ (pq + r) + 
m[y G (w - uq + vp) - z G (v + ur - wp)] = K f 

I y q + ( J x - ! z )P r - K (<7 r+ P ) + l yz (pq -r) + I xz (p 2 -r 2 )- 
m[x G ( w-uq + vp ) -z G (u- vr+wq)] = M f 

and 

I z r + ( I y- I x )pq-I xy (p 2 -q 2 )~I yz (pr + q) + I xz (qr- p) + 
m[x G (v + ur — wp) - y G (u — vr + wq)] - N f 

These preceding six equations are the most general form of Newton's laws for 
rigid body motion used today. They consist of a constant mass and inertia tensor, and are 
formulated in the body fixed frame. This set of equations can be simplified depending on 
the location of the local axes. If the axes are selected to coincide with the vehicle's 
principal axis of inertia, then the terms including the product of inertia become zero. 
This simplification is possible only if the xy-, xz-, and yz-planes are planes of symmetry. 
This is typically not the case. (Practically, most designs are only symmetric about the xz- 
plane. Further simplifications can be made, including locating the origin of the body 



(2.23) 

(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 
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fixed frame at the vehicle's center of gravity, however, these assumptions are not 
physically realizable.) 

The previous nonlinear system, equations 2.23-2.28, can be written in a compact 
vector form. Using the vector x the state vector of body fixed velocities, and defining 
Mrb as a rigid body mass matrix including translational and rotational inertial elements, 

m 0 0 0 mz G — my c 

0 m 0 ~mz G 0 mx G 

0 0 m my G -mx G 0 

M RE = ~ ITT'’ 

0 -mz G my G /„ I ^ I n 

mz G 0 -mx G I yx 1^ I yz j 

-my G mx G 0 l u I ^ I a 

and Crb(x) as a state dependent matrix containing the rigid body coriolis and centrifugal 



terms. 
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the equations of motion can be written in vector form as, 

M RB x + C RB (x)x = f. (2.29) 

The right hand term in Equation 2.29 is the vector of external forces and moments 
outlined in Table 2.1. These forces and moments come from hydrostatic and 
hydrodynamic forces due to gravity, radiation and excitation, viscous damping (lift and 
drag), and control inputs. The origin of these forces and their application to the 
developed equations of motion will be discussed in detail in following sections. 

D. RESTORING FORCES AND MOMENTS 

In hydrodynamic terminology, the gravitational and buoyant forces are called 
restoring forces. The weight, W, and buoyant, B, forces that act at the centers of gravity 
and buoyancy must be defined from static analyses. For submerged bodies the weight 
and buoyancy force vectors do not change with vehicle attitude. Assuming that weight 
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and buoyancy are fixed in relation to the body fixed frame, then the gravitational and 
buoyant forces may be expressed as, f w = 01 +0J + WK , and f B =01 +0J - BK . 



Since the weight and buoyancy terms in the applied forces act in the global vertical 
direction, they must be transformed into components in the vehicle fixed frame before 
they can be added into the equations of motion. Returning to Equation 2.4, it can be seen 
that the components acting along the vertical vehicle fixed frame are the third column of 
the transformation matrix. Therefore, the net vertical force components become. 



- sin 6 



/« 



= (W-B) 



cos 6 sin </> . 
cos 6 cos 0 



(2.30) 



The weight portion of the vertical force acts at the center of gravity of the vehicle, 
while the buoyancy portion of the vertical force acts at the center of buoyancy. Because 
these forces act in locations away from the body center they create a moment about the 



body center given by, 





- sin 6 




- sin # 


m g =Wp G x 


cos# sin <p 


-Bp B x 


cos # sin <p 




cos# cos (p 




cos# cos (p 



(2.31) 



This moment will be non-zero even if W and B are equal, since pc and pB are 
usually not collocated. For static stability it is necessary to locate the center of gravity 
below the center of buoyancy. The total vertical force vector can be written as, 

-(W-B)sme 
(IT -5) cos# sin (f) 

(W -B)cos0cos<p , . 

= -F(z), (2.32) 

(y G W - y B B ) cos 6 cos <p - (z c W - z b B ) cos 6 sin p 
— {x G W — x B B)cosdcos(p — (z G W — z s B)sin 6 
(x G W - x b B) cosd sm (f) + (y G W - y B B)sxn 6 

and is added to the right hand side of the equations of motion. 



fc = 



f* 



m 



8 _ 
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E. WAVE THEORY AND HYDRODYNAMICS 



1. Linear Wave Theory 



The simplest free surface wave formation, which nevertheless has great practical 
significance, is the plane progressive wave system where the water column is modeled as 
an inviscid, irrotational fluid in a gravity field. This motion is two dimensional, (x, z), 
sinusoidal in time with angular frequency (O, and propagates with a phase velocity c p such 
that to an observer moving with this speed the wave appears to be stationary. A 
Cartesian coordinate system is adopted, see Figure 2.2, with z = 0, the plane of the 
undisturbed free surface (still water level) and the z - axis positive upwards. The vertical 
elevation of any point on the free surface may be defined by a function z = 7](x,t). With 
these requirements, the free surface elevation must be of the general form 

Tj(x,t) = Acos(kx-0}t), (2.33) 

where the positive x-axis is chosen to coincide with the direction of wave propagation. 
Here, A is the wave amplitude, and the parameter 



k 




(2.34) 



is the wavenumber, the number of waves per unit distance along the x-axis. Clearly, the 
wavenumber can also be written as. 



k 




(2.35) 



where the wavelength L, is the distance between successive points on the wave with the 
same phase, see Figure 2.2. The solution of this problem is expressed in terms of a two 
dimensional velocity potential, which must satisfy Laplace's equation 

VV = 0, (2.36) 



and appropriate boundary conditions. Furthermore, the velocity potential, (p, must yield 
the wave elevation given by Equation 2.33 from 

77 = -—-^, at z = 0. (2.37) 

8 & 
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a Progressive Wave 
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Figure 2.2 Monochromatic Progressive Surface Gravity Wave [Rahman 1994] 

Equation 2.37 is the linearized dynamic boundary condition on the free surface and is an 
expression of the fact, through Bernoulli's equation, that the pressure on the free surface 
must be the same as the ambient atmospheric pressure. An appropriate boundary 
condition on the sea bottom is 

50 



5z 



= 0, at z = — H, 



(2.38) 



i.e., the bottom at depth H is a rigid impermeable plane. Finally, the free surface 
boundary condition is 



5 2 0 50 A A 

— = 0, at z-0. 

5 r dz 



(2.39) 



Equation 2.39 is a combined dynamic and kinematic surface boundary condition. 
The dynamic condition was discussed earlier, while the kinematic condition simply 
states, 

¥=¥■ < 2 - 4o) 

dt dz 

that the vertical velocities of the free surface and fluid particles are the same. Combining 
equations 2.37 and 2.40, Equation 2.39 is obtained, ignoring the small departures of the 
free surface 7] from the horizontal orientation of z = 0. 

From the requirements of the problem, it is clear that the velocity potential 0 must 
be sinusoidal in the same sense as Equation 2.30; therefore a solution of the form 
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<f>(x,z,t) = sm(kx-6X)F(z), (2.41) 

is attempted. Substituting Equation 2.41 into Equation 2.36, F(z) must satisfy the 
ordinary differential equation 

^^-k 2 F(z) = 0, (2.42) 

dr 

throughout the domain of the fluid. The solution to Equation 2.42 satisfying the bottom 
boundary condition is 

F(z) = Acosh(£(z + H )) . (2.43) 

Substitution of equations 2.41 and 2.43 into the surface boundary condition. Equation 
2.37, yields an important relationship between the wavenumber k and the frequency ox 

CO 2 = gk tanh(kH ) , (2.44) 



which is called the dispersion relationship. The surface elevation 7 ] follows from 
Equation 2.37 as, 

7j(x,t) = acos(kx-6X), (2.45) 

with the amplitude a, given by 



a=— cosh(kH). 

8 



(2.46) 



Substitution of equations 2.43 and 2.46 into the velocity potential function. Equation 
2.41, yields 



, . . ag cosh(£(z + H)) 

(p(x,z,t) = sm(fcx - ox) . 



(2.47) 



co cosh (kH) 

An underlying assumption associated with potential flow theory and Laplace's 
Equation, is that the velocity field can be expressed as the gradient of a velocity potential 
function. This allows the expressions for the velocities in the x and z direction to be 
given as 

dtp 
dx 
dtp 



u — - 



(2.48) 



w = - 



dz 



Using Equation 2.48, and the linearized form of the Bernoulli's equation. 
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dtp 

P=~P^ pgz. 



(2.49) 



the expressions for the fluid velocity and pressure fields are 



u — a 



gkcosh (k(z + H)) 



cos(kx - ax) 



at cosh (kH) 



gk sinh(£(z + //)) 
a> cosh (kH) 



sm(kx-ox) 



(2.50) 



cosh(k(z + H)) 

p = pga — — 

cosh (kH) 



cos (kx -ax)- pgz 



It can be seen from Equation 2.50 that the trajectories of the fluid particles are elliptical. 

There are several simplifications that may be made to the above-derived 
expressions for the cases of shallow (long waves) and deep (short waves) water. The 
shallow and deep water ranges correspond to H/L < tl TO and H/L > n, respectively. 
Over these ranges approximate expressions may be substituted for the hyperbolic 
functions that have been encountered. Table 2.2 summarizes these results. Figure 2.3 
depicts the comparison of water particle velocity between long and short waves. The key 
points that should be observed are that shallow water waves are non-dispersive, i.e., the 
vertical component of the wave particle velocity is linear in depth, and that the 
classification of shallow water depends on the ratio of water depth to wavelength. 

As an example, consider a one-meter high, monochromatic wave with a ten- 
second period (0.1 Hz), propagating in water six meters deep. The deep water and 
shallow water wavelengths are 156 meters and 76.7 meters, respectively. The ratios of 
water depth to wavelength are 0.038, for the deep water case, and 0.078, for the shallow 
water case. Referring to Table 2.2, it can be seen that neither the deep water nor shallow 
water simplifications may be used. In fact, the water depth must be reduced to 
2.45 meters before the shallow water equations are valid. Conversely, if the water depth 
remains at 6 meters, then the wave period must become greater than 15.6 seconds or the 
wave frequency less than 0.064 Hz. This indicates that only low frequency waves 
typically qualify as shallow water (non-dispersive) waves which becomes an important 
consideration in the proper modeling of the wave induced disturbances. 



26 





Deep water 
(short) waves, 
H/L> 1/2, 

kH>n 


Intermediate 
depth waves, 
1/20 <_H/L<A!2, 
71/ 10 <kH < 7t 


Shallow water (long) 
waves, 

H/L < 1/20, 
kH < tj/\0 


In all 
cases, 
9=kx-ca 








Profile, Tf 


77 = Acos# 


t\ = A cos# 


77 = Acos# 


Wave 

velocity, C 


C = — 

CO 


C=— tanh(«7) 

(O 


C = /gH 


Wavelength, 

L 


L 

In 


j gT 2 2 nH s 

L - 6 tanh( ) 

2 n L 


l=t^h 


Velocity 
potential, (/> 


Aco ~ 

<p = e sin# 

k 


Aco cosh(£(z + H )) . „ 

(p= sin# 

k cosh (kH) 


<t> = AgT sin# 

In 


Particle 
velocity, U 


u-Acoe cos# 


cosh (k(z + H)) 

u = Aco —cos # 

coshikH) 


Aco 

u = COS# 

kH 


Particle 
velocity, W 


iv = - A coe^ sin # 


A sinh (k(z + H)) . . 
cosh(fctf) 


w = — A co( 1 + — )sin # 

H 


Pressure, p 


p = pg(r]e kz -z) 


cosh(*(z + H)) 

P = pg(V z) 

cosh(kH) 


p=pg( n-z) 



Table 2.2 Summary Of Small Amplitude Linear Wave Equations 



As a side note, it should be pointed out that the wave height has no bearing in 
determining whether a wave is classified as long or short. The wave height is significant 
in order to determine the subsurface water particle velocities, and when the wave breaks. 

The plane progressive wave described so far is a single, discrete wave system, 
with a prescribed monochromatic component of frequency co and wavenumber k, moving 
in the positive ^-direction. More general wave motions, which are not monochromatic. 
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(a) Shallow water 




( b ) intermediate depth 




( c) Deep wa ter 



Figure 2.3 Particle Orbits And Variation Of Particle Velocity Amplitude With Depth 

[Sarpkaya 1981] 



can be obtained by superimposing plane waves of different frequencies and 
wavenumbers. The elevation of the sea surface ij(t) can thus be described as the 
superposition of an infinite number of sinusoids of the form: 



*7(0 = 2X c °s(K x ~ W + <P n ) =2X • 



(2.51) 



n = 1 



n = 1 



The equations for the horizontal and vertical velocities, and the dynamic pressure are, 

co n coshk n (H + z) 



u(t) = £ 

n = 1 

w(r) = t 



n=\ 



sinh k„H 



co n sinh k n {H + z) 
sinh k„H 






7 n ; 



Pit) = PsYj 



n-\ 



a> n cosh k n (H + z) 
cosh kH 



In - P8Z ; 



(2.52) 

(2.53) 

(2.54) 
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respectively, where the parameters/variables in equations 2.51-2.54 are identical to those 
in the monochromatic case. 

2. Hydrodynamics 

[Newman 1977] has shown that the total velocity potential </>, may be written as 
<t> = <t> w +<P D +<Pr, the linear sum of three components. These three components come 
from wave, diffraction and radiation potentials where, 

• <pw is the incident regular wave velocity potential; 

• <pD is the diffraction potential caused by reflection when the vehicle is considered 
restrained from motion; and 

• (j>R is the radiation potentials in 6 DOF caused by forcing the vehicle to oscillate with 
wave excitation frequency, when there are no incident waves present. 

Linear wave theory implies that the wave induced forces and moments acting on a 
vehicle come from the superposition of the radiation induced forces and moments and the 
excitation forces and moments. 

Radiation induced forces and moments act on the vehicle when the vehicle is 
forced to oscillate with the wave excitation frequency without incident waves present. 
The forces due to wave radiation are classified as restoring, added mass and potential 
damping forces. 

Excitation forces and moments act on the vehicle when the vehicle is restrained 
from motion and incident waves are present. These forces and moments caused by wave 
excitation are classified as Froude-Kriloff ( FK) and diffraction forces. The FK forces and 
moments are found by integrating the pressure distribution over the vehicle cause by the 
undisturbed wave field, while the diffraction forces and moments are determined by the 
pressure distribution created when the waves are reflected from the vehicle. 

As a general summary, the unforced equations of motion for any stationary body 
in waves, can be given as 

(M RB +A( 0 )))fj + B(co)f}+C7]-(M FK + \{co))rj f -B(a>)fj f -Cjj f =0, (2.55) 

where tj is a vector of linear and angular displacements. Stationary in this sense implies 
that the vehicle is non-maneuvering and that its only motion is that caused by the body’s 
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interaction with the waves. As seen, the added mass and damping matrices, A( a) and 
B( co) respectively, are functions of the incident wave frequency. However, if the vehicle 
is moving with some forward velocity U, the wave frequency, ox is not the frequency 
encountered by the vehicle. The encounter frequency, ox, can be expressed as 



In this equation, which is based on the Doppler effect, /? is the heading angle between the 
vehicle and the wave propagation direction and k is the wavenumber. More detailed 
information on the Doppler effect maybe found in Appendix B. 

As a simplification, if the body is small in comparison to the wavelength, it is 
totally submerged and neutrally buoyant with homogeneous mass distribution, then 
Equation 2.55 can be simplified to yield 



where f] r = fj — TJ f , is the relative velocity of the fluid over the vehicle. It is necessary to 

point out that the above equation was based on the assumption of inviscid flow. Due to 
this fact, viscous damping terms (skin friction and drag), must be taken into account for a 
vehicle operating in a real fluid, to complete the model. 

In 1950, Morison developed an expression for the horizontal force on a 
cylindrical pile subjected to waves. His work showed that the elemental force dF on a 
vertical strip of a cylinder dz, may be written as 



In this expression, coined “ Morison' s equation ", D is the characteristic diameter, u f 

and uj are the horizontal component of the undisturbed fluid acceleration and velocity, 
respectively, and Cm and Cd are mass and drag coefficients which may be determined 
experimentally. This equation has two parts, the first term representing an inertia force 
proportional to the accelerating fluid acting on the pile (a mass force), and the second 
term, a nonlinear drag term proportional to the sign squared fluid velocity (a drag force). 

In practice, Morison's equation can only be applied to small volume bodies. By 
small volume bodies, it is meant that the characteristic cross sectional dimension of the 



0) e =Q)+kU cos/?. 



(2.56) 



(M KB +A(Q} e ))ii r +B((o e )7i r +C7i r = 0, 



(2.57) 




(2.58) 
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body is small relative to the wavelength [Morison 1950], For vertical piles, small volume 
applies only if L>5D, where D is the pile diameter, refer to Figure 2 4 As shown, for 
small vehicles operating in shallow water, inertia and drag forces are dominant, where as 
reflection and diffraction effects can be considered unimportant 

It is known that the ratios of wavelength and wave height to the characteristic 
diameter are key parameters in predicting the load regime of waves acting on a structure 
[Faltinsen 1990] These regimes are also depicted in Figure 2,4. 

For a stationary object in a simple harmonic flow, the total force can be expressed 



as 



F T = F d sin tut sin (Ot\ f F, cos cot , 



(2.59) 



where Fd and F/ represent the maxima of the drag and inertia force components, 
respectively. [Dean 1984] divided the flow force regimes into two areas; one where the 
inertially derived component dominated the total force, and one where both drag and 
inertial effects are important. This expression is represented as 




Fn 



F d + 




2 F d <F,' 
2F d > F, “• 



(2.60) 




Figure 2.4 Load Regimes 



31 



The significance of the above expression is that the maximum force on the body is not 
affected by additional drag force until the amplitude of the drag is at least one half that of 
the inertia force. For oscillating flows, such as those caused by waves, while even small 
amounts of drag may be important when considering the shape of the load function on a 
stationary body, the peak amplitude of the force is only affected when the drag 
component is greater than one-half the inertia force. Extensive experimental verification 
of Morison’s approach to force modeling, and the evaluation of the frequency dependent 
nature of the drag and added mass coefficients was given by [Sarpkaya 1975]. 

When the wave field causes motion of the water particles at the vehicle's 
operating depth to be of the order of one diameter or less, it is expected that the 
predominant hydrodynamic force on the AUV due to the wave disturbance would be 
inertial in nature. Since the AUV operates below the surface, it is not the wave height to 
vehicle diameter ratio that is of concern, but more appropriately the double amplitude of 
the water particle motion at the vehicle operating depth compared to the vehicle 
characteristic length of interest. While this analogy is approximate in nature, it does 
provide a means to predict which hydrodynamic forces may be of concern when 
estimating the total load on a vehicle. These above concepts may be used to estimate the 
dominant forces acting on an underwater vehicle subject to waves, and therefore assist in 
the sizing of the propulsion system. 

F. HYDRODYNAMIC FORCES AND MOMENTS 

Hydrodynamic forces and moments are the result of body/fluid interactions. The 
forces and moments on the body arise from the modification to the pressure distribution 
summed around the surface area of the body. This modification to the pressure field can 
only arise from relative velocity and acceleration between the body and fluid. Therefore, 
for the purposes of this discussion it is necessary to re-define the body fixed velocity 
vector x, in terms of a relative body fixed velocity vector x r , where 
x r = [u r , v r , w r , p, q, r ] . Also at this time it is convenient to define a globally based fluid 
velocity vector Uf, where t//= [U/, Vf, Wf, 0, 0, 0] r . Since it is assumed that the fluid 
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velocity is irrotational, no changes to the angular rate terms are necessary in the body 
fixed velocity vector, and no angular rates are present in the fluid velocity vector. 

1. Radiation Induced Forces and Moments 

a) Added Mass 

Like the rigid body kinematics, it is desirable to separate the added mass 
terms into terms which belong to an added mass matrix Mam and a matrix of Coriolis and 
centrifugal terms Cam(x). For underwater vehicles this implies that the added mass 
forces and moments can be written as: 

f am AM x r —C AM (x r )x r (2.61) 

where f A M=[ X a,Y a ,Z a , K&M A , N a] -[/a, m A ] is the total added mass force and moment 
vector. 

The added mass terms represent the inertial reaction of fluid particles 
surrounding the submerged body that are accelerated with it. Any motion of the vehicle 
induces a motion in the otherwise stationary fluid. In order to allow the vehicle to pass 
through the fluid, the fluid must move aside and then close behind the vehicle. As a 
consequence, the fluid passage possesses kinetic energy that it would lack if the vehicle 
was not in motion. 

[Lamb 1932] gives the following expression for the fluid kinetic energy, 
Ek, which may be expressed in a quadratic form of the body axis velocity vector 
components; 

E*=-{pl<t’^f; = \x, T M AM x r . (2.61) 

In Equation 2.61, Mam is a 6x6 added mass matrix. For a rigid body moving in an ideal 

T 

fluid the added mass matrix is symmetrical, i.e.. Mam = Mam ■ In a real fluid these 36 
elements may all be distinct. [Wendel 1956] has shown that the numerical values of 
added mass in a real fluid are usually in good agreement with those obtained from ideal 

T 

theory. For a body possessing vertical plane symmetry only, and applying M a m = Mam , 
the added mass matrix is written as 
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X, 
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X, 


0 
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r, 


0 




** 
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2, 


0 




0 


0 


Yp 


0 




0 


Kr 


** 


0 




0 




0 


0 


Y r 


0 


K r 


0 





The notation of [SNAME 1950] is used in this expression. This notation indicates the 
degree of freedom on which the hydrodynamic added mass force acts, as well as the 
cause of the force. As an example, Y^u. is a force acting along the body fixed y-axis due 



to an acceleration u in the x-direction, and can be thought of mathematically as 



Y i = — . This definition implies that the hydrodynamic derivatives corresponding to the 
du 

diagonal of the added mass matrix will all be negative. 

The added mass terms are obtained from potential theory. This theory 
assumes an inviscid fluid, no circulation and that the body is completed submerged in an 
unbounded fluid. The last assumption is violated at the seabed, near underwater objects 
and at the surface. [Milne-Thomson 1968] has shown that the expressions for the force 
and moments resulting from added mass effects can be found by applying Kirchhoffs 
equations. 



d_ 

dt 

d_ 
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du 



r dE t 



+ fi)X- 



d(0 



+ C0X- 



BE, 

du 
dl ? 



= -/„ 



- + UX- 



dl u 



(2.63) 



= -m 



d(0 du 

in vector form. Expanding Equation 2.63, the expression for the added mass terms 
associated with the x-direction is 



X A = X-u + X^w+uq) + X-q + Z.,wq + Z-q 2 

-Y.vr-Y^rp-Y.r 2 



(2.64) 



(Derivation of the added mass terms associated with the other degrees of freedom is left 
to the reader.) Many of the added mass derivatives contained, in the general form are 
either zero or mutually related to another term when the body has various symmetries. A 
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more detailed discussion of added mass terms applied to an underwater vehicle, is found 
in [Humphreys 1978]. 

Extracting the added mass derivatives corresponding to the velocity 
coupling terms from 2.63 yields; 



where, 
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b ) Wave Radiation or Potential Damping 

The contribution from the potential damping terms compared to other 
dissipative terms like viscous damping terms are usually negligible for underwater 
vehicles operating at great depth. Nevertheless, underwater vehicles operating in shallow 
water near to the free surface should consider potential damping effects, especially, those 
underwater vehicles that tend to have a non-streamlined body, i.e., vehicles build with 
sensor and equipment mounted in such a manner that the equipment causes the vehicle 
not be streamlined. The linear potential damping can be modeled as 

f d =~D d x r . (2.65) 

The linear damping matrix D d is a positive definite matrix of linear damping coefficients. 
These linear damping terms are small when compared to the viscous forces, and therefore 
are often included in the viscous drag forces. 
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2 . 



Excitation Forces 



When applying potential theory, the fluid motion was assumed to be irrotational. 
This implies that only the linear velocity components, V /, of the fluid are considered 
when determining the excitation forces. In linear theory, the wave induced forces and 
moments acting on a vehicle are considered to be the sum of the radiation induced forces 
and moments and the excitation forces and moments. For nonlinear theory, this is not the 
case. The forces and moments due to radiation and diffraction are nonlinear functions of 
the relative velocity and acceleration between the vehicle and the fluid. 

a) Froude-Kriloff Forces 

The Froude-Kriloff force and moment vector can be expressed as 

f FK = M FK^f 5 ( 2 . 66 ) 

where M F k can be interpreted as the Froude-Kriloff ( FK) mass matrix and u f is the fluid 

acceleration expressed in body fixed coordinates. Coriolis and centrifugal terms will not 
appear in the general expression for the FK forces and moments since it was assumed that 
the rotational fluid motion was zero. This mass matrix can be determined by computing 
the mass of the fluid displaced by the vehicle and substituting this value into the rigid 
body mass matrix in place of the vehicle mass. Also, it should be recognized that this 
excitation force acts at the center of buoyancy not the center of gravity. The mass and the 
moments and products of inertia for the FK mass matrix are 

m = pV =B/g, T„ = 2" dm, {y 2 + r ) and 1 ^ = I yx = (xy ) . 

i=i 

Substituting these values into the rigid body mass matrix yields the FK mass matrix, for a 
small volume completely submerged body. 
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[Sarpkaya 1981] and [Newman 1977] have shown that the body fixed 
inertia force to which a symmetric body moving in an unsteady flow field is subject, may 
be written as 



/ = />V(/+C,) 



du. du. 



-pVC A u 



(2.67) 



In this equation, it is important to note the presence of a buoyancy-like force that is 
proportional to the displaced volume of the body as well as the presence of some 
convective terms . The terms that include the displaced volume, when grouped together, 
represent the added mass and Froude-Kriloff forces. The convective terms represent the 
forces from the spacial changes of the unsteady flow field over the body, and for large 
bodies, may be significant. [Silvestre 1998b] has shown that by comparing the force 
contribution of the convective terms to the total inertia force exerted on a rigid body 
subject to wave disturbances, that the forces due to the convective terms are small and for 
all practical purposes may be neglected. 

3. Viscous Damping Forces and Moments 



a) Drag Force 

The drag effects on an underwater vehicles are mainly caused by 

• Linear skin friction due to laminar boundary layers; 

• Quadratic skin friction due to turbulent boundary layers; and 

• Quadratic drag due to vortex shedding (Morison's equation). 

The viscous damping forces and moments will be functions of the relative fluid motion. 
In the range of Reynolds numbers in which underwater vehicles typically operate, flow is 
turbulent, therefore the drag force is approximated by the square law resistance arising 
from Morison's equation. Referring to Equation 2.55, the quadratic drag force in the 
x-direction can be expressed as 

/ =^P C D Au r\ u r\ ( 2 - 68 ) 

where A is the projected cross-sectional area, Co is the drag-coefficient based on the 
representative area, p is the fluid density and u r is the relative longitudinal velocity. 
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A generalization of Morison's equation could be to use a truncated second 
order Taylor series expansion to describe the viscous damping in 6 DOF. This suggests 
that the viscous damping forces could be written as 

/„ = ~Dj x r -D q (x r )~ higer order terms . (2.69) 

This expression is a combination of linear and quadratic drag components. The linear 
portion D,x r , can be modeled as a positive definite matrix of linear damping coefficients 

multiplied by the body fixed relative velocity vector, while the quadratic portion D q (*,), 

is a nonlinear matrix incorporating the contributions due to skin drag and vortex 
shedding. 

It is quite complicated to determine the hydrodynamic coefficients 
associated with the quadratic portion, especially for high angles of attack. However, a 
simple, but fairly accurate method of modeling these viscous drag forces is to include the 
linear drag components associated with the diagonal of D t x r , and to use a cross flow 

integral to account for the D q (x r ) terms. The general cross flow drag expression for a 
body of revolution in the y-direction (sway) is given by 

Y drag--\p \ [C Dy D(x)(v + xr) 2 + e Dz D(x)(w -xq) 2 ] dx . (2.70) 

In this expression, the minus sign is present because the drag force opposes the motion, 
and the cross flow velocity is given by 

U cf (x) = [(v + xr) 2 +( W -xq) 2 ] u2 . (2.71) 

Equation 2.70 is valid for a vehicle that has a hull form consistent with a 
body of revolution, however, for underwater vehicle's with different shapes this equation 
must be modified. As an example, consider the box shape hull form of the NPS Phoenix 
AUV [Marco 1996]. The cross flow drag expression in the sway direction for this shape 
vehicle becomes 

Y drag =-~P \ C D r h (x)(v + xr)|(v + xr)\dx . (2.72) 

2 - 1/2 
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As depicted in Figure 2.5, this expression reflects the drag force over the entire range of 
angles of attack. 

b) Lift Forces 

Lift is a force acting on a body in a direction perpendicular to that of the 
flow of the fluid. As relative motion is created between the underwater vehicle and its 
fluid environment, the vehicle body experiences lift forces similar to those experienced 
by an airfoil. The expression for the lift generated by an airfoil is given as 



where A is the projected cross-sectional area, Cl is the lift-coefficient, a is the local angle 
of attack, p is the fluid density, and u r and v r are the relative velocity components in the 
body fixed x and y directions. See [Hoemer 1975] for further detail. 

The lift forces can be modeled as a state dependent matrix multiplied by 
the respective relative velocity and is given by 



The parameters/coefficients used in this matrix, like the drag coefficients are difficult to 
predict and vary with the shape of the body and the angle of attack of the fluid. As with 
the drag forces, determination of these constant coefficients is typically done using a 3-D 
CFD solver, with the constant coefficients validated over a range of angles of attack. 

G. CONTROL FORCES 

Small underwater vehicles are usually maneuvered with thrusters and control 
surfaces. Thrusters are effective only at low forward vehicle speed due to the fact that 
the action of the thrust is a nonlinear function of the relative velocity of the vehicle. 
However, control surfaces are effectively used in maintaining heading as well as trim and 
depth changing maneuvers. The reason for this is that the force generating capacity of 
the control surfaces is dependent on the speed of the vehicle (the lift force is proportional 
to the square of the velocity). The control forces and moments can be described by 




(2.73) 



f L =-D L (x r )x r . 



(2.74) 



fc =B(x r ,u comrol )u 



control 



control 



(2.75) 
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Figure 2.5 Normalized Drag Force vs. Angle of Attack 



where u CO ntroi is an input vector and B(x n u co mroi) is a state dependent input matrix. This 
input matrix contains the necessary coefficients to model the forces developed by the 
control actuators. The total force/moment vector fc, is the sum of the propulsion force 
vector,^,, due to thrusters and propellers, and the actuator force vector, /& due to fin and 
rudder deflection. 

1. Propulsion Forces 

The derivation of a steady-state hydrodynamic model for propellers operating in 
an incompressible fluid can be found in many introductory fluid texts see [Lewis 1988, 
Newman 1977 and White 1986] for examples. These models are based on large open 
propellers. The use of small thrusters for control of underwater vehicles is an area of 
current research. This is because the vehicles are small, require fast response and are 
required to conduct dynamic positioning maneuvers. [Yoeger 1991] developed a lumped 



40 



parameter model that improved on the popular notion that for a given unit the thrust and 
input torque are related to the square of the propeller rotational rate and the angle of 
advance. He also showed that by accounting for thruster dynamics, improvements in 
position control could be obtained. Yoeger's work introduced the concept that fluid 
momentum considerations resulted in a time lag in the thrust response to a step input. 
Although his work improved the modeling of thrusters, it left room for significant 
improvements. 

[Healey et al., 1995], improved on Yoeger's work by providing a generic thruster 
model that considered propeller thrust and torque as a mapping linked to lift/drag force 
variations caused by changes in the local angle of attack of the propeller blade. They also 
were able to associate the lags and overshoots in the thrust response to lags in the 
development of the local angle of attack and dynamic development of the blade pressure 
distributions. Research by [Whitcomb 1995] and [Bachmayer 1998] has provided further 
experimental validation that the four-quadrant model proposed by Healey is valid for 
small ducted thrusters. 

The model proposed by the researchers in the above paragraph can be best 
expressed as a first order differential equation of the form 




where F is the propulsive force imparted to the vehicle. The parameter r is the time 
constant associated with the force lag, the second term containing y is a thrust reduction 
cause by the change in local angle of attack as the propeller advances through the water 
and the last term containing (3, is the thrust to rotation rate mapping. In standard 
propeller design terms, the parameter yean be related to J, while /? can be related to K T . 

Referring to Figure 2.6, it can be seen that the thrust coefficient Kj is a function of 
the propeller speed of advance J. The value of the non-dimensional thrust coefficient at a 
zero speed of advance is associated with the bollard pull condition. This relates the thrust 
to the rotational rate of the propeller as 

T = K t pD A n\n\, (2.77) 
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where D is the diameter of the propeller and n is the rotational rate. 

For any other operating condition, it can be seen that the thrust coefficient is 
reduced as the propeller moves through the fluid with some speed of advance. The thrust 
for this operating condition is given by 

T = K T pD*n\n\ + y o JpD 4 n\n \ , (2.78) 



where y 0 is the slope of the Kj curve at the particular operating condition of interest 
(negative over the entire curve), and J is the non-dimensional speed of advance. 
Substituting the expression for J, shown in Figure 2.6, the parameters (3 and y from 
Equation 2.76 can be related to the parameters in Equation 2.78 as 



t 6 = K T pD 4 

r=Y 0 pD 3 



(2.79) 




a 

sc 



Figure 2.6 Thrust and Torque Coefficients versus Angle of Attack [Lewis 1988] 

The force/moment vector caused by the propulsion forces due to the set of 
thrusters and main propellers for the NPS Phoenix AUV is, 
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(2.80) 



0 





where the x and y distances are all positive values measured from their respective plane 
of origin. 

2. Actuator Forces 

Control of a small AUV at anything other than slow speed must be accomplished 
by using control surfaces since the effect of non-propulsion thrusters decreases as the 
forward speed of the vehicle increases. These control surfaces are comprised of 
fins/planes and rudders. The forces and moments exerted by these actuators are derived 
from airfoil theory and are composed of lift and drag components. With the exception of 
the longitudinal direction, the force/moment applied to the vehicle is directly proportional 
to the amount of angular deflection of the control surface. In the longitudinal direction 
the force exerted by a deflection of a surface amounts to an additional drag force on the 
vehicle. The vector of forces and moments caused by these surfaces is 



The modeling represented in Equation 2.81 is for a vehicle equipped with a 
standard fin arrangement. By this it is meant that each pair of control surfaces, stem 
planes, bow planes or rudders, move together and are not independently controlled. If the 
vehicle is equipped with either an X-brace configuration[Humphreys 1994], or the 
actuators in a standard configuration are allowed to be moved independently, then the 




(Y Ss S sr +Y Sb S br )u\u\ 
(Z Ssp S sp +Z K S bp )u\u\ 



. (2.81) 



0 



( M S sp S sp +M sJ bp >\ U | 
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modeling of these control forces will change. Specifically, there will be the ability to 
actively control roll and there will be some level of redundancy for each of the other 
control modes. For a detailed description of the origin of these control forces, the reader 
is referred to [Hoemer 1965, Hoemer 1975 and Lewis 1988]. 



H. 6DOF EQUATIONS OF MOTION 
1. No Fluid Motion 



Using Equation 2.29, and substituting the external forces discussed in the 
previous sections, the 6DOF EOM can be written as 

M rrX +C RB (x)X = f G + /am + fd "b f FK + fv ■b f L "b fc (2 82) 

Expanding Equation 2.82 by substituting the appropriate expressions for the 
force/moment vectors, and recognizing that without fluid motion x = xr and uf = 0, the 
following equation is obtained 

Af rrX + C Rg(x)X + F (z) + M am % + ^ AnfX^X + 

D d x + D,x + D q (x) + D L (x)x = B(x, u con[rol 

control 



(2.83) 



As seen. Equation 2.83 contains a mixture of coordinates; both body fixed and 
global. To use a system of equations expressed in this form for simulation studies or 
control system development, this system must be augmented with the relationships 
between the various coordinate frames. The link between the global and body fixed 
coordinates is accomplished by augmenting Equation 2.83 with equations 2.8 and 2.11. 
This system is represented by 

[Af/fB b A 1 AM ]x + [C ^(x) + C Am (x)]x + 

[D d +D,+ D l (x)]x + D q (x) + F(z) = B(x,u comrol )u conlrol (2.84) 

z = g(x,z) 



2. Modifications To Account For Fluid Motion 



In the case where fluid motion is present. Equation 2.83 is represented as 
A t rrX + C RB (x)x + F (z) + M AM x r +C A Jx r )x r +D d x r 

~ A^ fk^ f "b F) \X r + Dq (x r ) + D l ( x r )x r =B(x r ' H control control 
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This equation, as was Equation 2.83, is a mixture of various coordinate frame 
variables; body fixed, body fixed relative, global and fluid. To solve this system of 
equations, the equations must be expressed in variables that can be related to each other. 
Since, for this case, the vehicle is considered to be in an unsteady fluid referenced to the 
global frame (refer to Figure 2.1), the logical choice for variables is body fixed relative 
xr, and global, z. 

Remembering that relative velocity is defined as Xr = x - u h and that the body 
fixed fluid velocity can be expressed as 



u f = 



T t 0 
0 7\ 



U f =TU f 



( 2 . 86 ) 



Equation 2.85 can be modified so that it is expressed in body fixed relative and global 
variables. Manipulating Equation 2.85 results in 

[M rb +M am ]x r +[C RB (x r )+C Am (x r )]x r +[D d +D v (x r ) + D L (x r )]x r +F(z) 



= —M 



RB 



d(TU f ) d(TU f ) 

' -C m (u f )TU f -—L. + B(x r ,u M )u ml 



dt 



dt 



(2.87) 

Looking at Equation 2.87 and recalling that the fluid was defined as irrotational, it can be 
seen that the Crb(u/) term on the right hand side of the equation is zero. The other two 
terms, on the right hand side, containing the rigid body mass matrix Mrb, and the 
Froude-Kriloff mass matrix Mfk, are excitation forces resulting from the fluid motion. 

As seen in Equation 2.87, the Froude-Kriloff excitation forces and moments are 
functions of the weight and buoyancy mismatch (Vf-5), the fluid velocities and fluid 
accelerations, expressed in body fixed values. For a neutrally buoyant vehicle, where 
W=B, the wave excitation forces do not present themselves in the translational equations 
of motion, however, they still provide excitation moments to the rotational equations. 
The reason behind this is because the fluid components (acceleration and velocity) act at 
the vehicle center of buoyancy, while the body inertial acceleration reaction force acts at 
the vehicle’s center of mass. 

At this point the 6DOF EOM representing the vehicle dynamics has been 
modified to account for a moving fluid by representing the body fixed velocities in 
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relative terms. Equation 2.87. However, as with the system of equations for the case of 
no fluid motion, the system in Equation 2.87 must be augmented to provide the necessary 
link between the global and body fixed velocities. In order to account for the fluid 
motion, either wave induced or steady current, Equation 2.8 is modified and is 
represented as 



~x 




~u r ~ 




U f 


Y 


= T;\<t>,6,yf) 


v r 


+ 


Vf 


Z 




_ W r_ 




Wf 



( 2 . 88 ) 



where the fluid velocity is represented in the global frame. Combining equations 2.87, 
2.88 and 2.11 into a system of equations results in the necessary equations to describe the 
motion of a small underwater vehicle subject to shallow water waves. 

It was shown earlier in this chapter that the total time derivative is composed of a 
rate of growth term and a rate of transport term. In the case of the time derivative of the 
rotation transformation matrix T, since there is no translation, the time derivative can be 
expressed solely as axT . [Healey 1992a] and [Fossen 1994] have shown that the 
cross-product operation (OX), can be represented as a matrix by defining the screw 
symmetric matrix, 



S(6>) = 



0 — r 

r 0 

-q p 



<7 

~P ■ 
0 



(2.89) 



The elements of S( to) are based on the consideration that both vector-matrix quantities in 
the cross-product operation are represented in the same coordinate frame. In the case 
where T is the transformation matrix from global to body fixed coordinates, the 
expression for the rotation matrix time derivative is 



S t (0))T, 0 

0 T 2 



(2.90) 



S T (a>) is the transpose of the matrix represented in Equation 2.89 since T is not in the 
same coordinate frame as ax Using Equation 2.90, and the fact that the coriolis matrix, 
Crb(u /), is null, the vehicle dynamics equations expressed in matrix form are 
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(2.91) 



The complete expanded 6DOF EOM, including physical parameters for the NPS Phoenix 
AUV, are outlined in Appendix A. 

I. DEVELOPMENT OF LONGITUDINAL SURGE MODEL 

Restricting the motion of the vehicle to surge only, the significant 
motions/quantities that must be incorporated to effectively model the vehicle in the 
longitudinal direction are, the surge velocity u r , and the global position X. This 
restriction simplifies the twelve previously developed equations to a system of two 
non-linear equations of motion. Based on the NPS PHOENIX AUV equations of motion, 
see Appendix A, this reduced set which models longitudinal surge motion is, 



This set of equations, with a slight modification, is also valid for modeling relative 
motion. By representing all variables in body fixed quantities, the set of equations for 
relative positioning is, 



where the position x, is measured relative to the vehicle fixed frame. 

The purpose of this longitudinal surge model is to allow for the development of a 
surge controller that will allow a vehicle to hold position in the presence of waves. To 
complete this model, the propulsion system dynamics must be included. Therefore, the 
system given in Equation 2.93 must be augmented with Equation 2.76. Augmenting 
Equation 2.93 with the propeller force equation, and assuming a neutrally buoyant 
vehicle, the set of equations to be used as the basis of the longitudinal surge controller 
becomes, 




(2.92) 



X =u r +U f 



(m-X.)u r +X uH u r \u r \ + 




(2.93) 



x = u r +u f 
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where the parameters a, /?, y and r must be determined through system identification. 
The process by which these parameters are determined will be presented in Chapter IV. 
Note, it should be pointed out, that in the form used in Equation 2.94, F becomes a 
generalized force with units the same as an acceleration rather than a direct thrust value. 

J. SUMMARY 

This chapter has presented the kinematic and dynamic relationships used in 
modeling a small underwater vehicle operating in a shallow water wave environment. It 
discussed the various external forces and moments that act on an AUV. The chapter 
concluded with the development of the one degree of freedom (1DOF) surge model that 
will be used as basis for a controller that will allow an AUV to station-keep in the 
presence of waves. 



IU. DISTURBANCE ANALYSIS 



A. INTRODUCTION 

In general, the disturbances that act on an underwater vehicle can be placed in 
three categories and described as follows: 

• Additive disturbances are external forces and moments which act additively 
on the vehicle. By including their effects, the total description of the vehicle 
model is extended by additional states (e.g., current, waves, and wind). 

• Multiplicative disturbances affect the dynamics of the system (e.g., the depth 
of the water, load conditions, trim, and speed changes). These disturbances 
can be regarded as time variant. 

• Measurement disturbances are due to incorrect measurements (e.g., noise on 
the vehicle's sensors). 

For this dissertation, only the additive disturbances will be taken into account. 
From the class of disturbances causing additive effects on an underwater vehicle in 
shallow water, only waves and current will be considered since they are most dominant. 
The ability to control a vehicle is known to be significantly affected by its environment. 
Since the modeling of external disturbances, especially waves, is rather complicated, 
many attempts to design control systems have suffered. 

A general assumption that is used in the modeling of AUVs is that forces and 
moments, which are added to the “calm sea ’ 4 model, can model environmental induced 
disturbances. This procedure was outlined in the previous chapter. This method, using 
the principle of superposition, is a good approximation for most marine control 
applications; however, it should be noted that for large general maneuvers it is not 
expected to be valid. 

This chapter will begin with a discussion of the stochastic nature of sea waves, 
including a description of several empirical relationships that can be used to represent the 
spectral content of a wave field. This will be followed by an overview of state space 
representations and recursive modeling of wave disturbances with specific application to 
control design. Finally, a methodology for the use of empirically derived spectral 
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relationships as well as measured wave elevation time series in distributed simulations 
will be outlined. 

B. STATISTICAL DESCRIPTION OF SEA WAVES 

The proper characterization of the real sea surface is difficult to obtain. With 
respect to the design of a control system, it is sufficient to assume a simplified description 
of the sea by considering only unidirectional linear waves. Based on this assumption, and 
on the superposition principle, simplified models can be determined. Concerning the 
modeling of waves, there exist two typical approaches: regular waves and irregular 
waves. 

Regular waves can be represented as a simple two-dimensional, sinusoidal wave 
train over an infinite water surface with infinite depth. This interpretation is based on the 
empirical observation, that the motions generated by waves have strong periodic 
components, see [Kallstrom 1979] for details. The characterization of the sea level as a 
train of regular waves, is however, an approximation, which is not necessarily accurate. 

Irregular waves allow the stochastic nature of waves to be taken into account. 
According to this method, the sea level can be modeled as either a superposition of a 
large number of regular waves of different amplitudes, frequencies, phase angles and 
directions of propagation, or as a narrow band stochastic process. In the case where the 
sea level is regarded as a stochastic process, the spectral and probability density 
description should be available, thus the model for the variation of the wave elevation can 
be determined from its spectral density function. The irregular waves are considered 
probabilistically with respect to amplitude and wavelength. Since the origin of waves is 
usually due mainly to the wind, the frequency and the steady state amplitude of the waves 
depend on the mean value of the wind speed. In this and the following sections, the 
- stochastic characteristics of waves are considered. 

The wave elevation of a long-crested irregular sea propagating in the positive 
x- direct ion can be written as the sum of a large number of wave components represented 

by. 
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(3.1) 



rj(x, 0 = X A > sin (M “ co t t + <f > i ) 

i=i 

where. At is the wave amplitude, cot is the wave frequency, (p t is a random phase angle and 

2k 

*' - r (3.2) 



is the wave number, with L,- as the wavelength. The wave amplitude can be expressed as 
a wave spectrum S( co), by 

A 2 = 2S(a> i )A.ct> (3.3) 



where, Ary is a constant difference between successive frequencies. The instantaneous 
wave elevation has a Gaussian distribution with a zero-mean and a variance a 2 defined 
as 



a 2 =\~S(co)do). (3.4) 

1. Wave Spectral Densities 

Several formulations of wave spectral densities have been proposed. The four 
spectra commonly encountered in practice are: 

• The Amplitude Spectrum. In this spectrum, the ordinates of the spectral 
density are proportional to the amplitude squared of the component waves. 

The area under the spectrum curve, Sa , is proportional to twice the average 
energy of the record. 

• The Energy Spectrum. The ordinates of the spectral density are proportional 
to half the amplitude squared of the component waves. The area under the 
spectrum curve, Se, is proportional to the average energy of the record and 
equals SJ 2. 

• The Height Spectrum. The ordinates of the spectral density are proportional 
to the height squared of the component waves. The area under the spectrum 
curve, Sh, equals 4S A . 

• The Double Height Spectrum. The ordinates of the spectral density are 
proportional to twice the height squared of the component waves. The area 
under the spectrum curve, S 2 H, is therefore equal to 85,4. 

Graphical representations of these spectra are shown schematically in Figure 3.1. 
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2. Statistical Description of Wave Amplitudes 

The probability density function of wave amplitudes with a narrow band spectrum 
can be expressed by the Rayleigh distribution, 

P(r) = ^-e~ r2,?2 (3.5) 

r 

where p(r)dr is the probability that a wave amplitude (r) lies between r and r+dr, and r 2 
is the mean square value of the wave amplitudes in the record, [Longuet-Higgins 1953]. 
It has been shown through the use of histograms that actual wave amplitudes closely 
follow this theoretical distribution, therefore, with the use of Equation 3.5, some 
quantitative statistical results may be formed. 
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Figure 3.1 Wave Spectral Density Comparisons [Berteaux 1976] 



a) Most Probable Wave Amplitude 

The most probable wave amplitude is the value of r for which 

~r p( r )~ 0 • (3.6) 

dr 

Differentiating the probability density function. Equation 3.5, and equating it to zero as 
indicated in Equation 3.6, yields an expression for the most probable wave amplitude r m , 
as 

r m = 0.707 (3.7) 
where is the root mean square value of the wave amplitudes in the record. 
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b) Mean Amplitude 

The fraction /(0 </< 1) of wave amplitudes larger than a given amplitude 
r Q is represented by 

f=r P (r)dr. (3.8) 

ir o 

The average of the / highest amplitudes can be found from the integration 

r f =^-\~rp(r)dr. (3.9) 

/ Jr ‘ 

The mean amplitude of all the waves in the record is then obtained when/= 1 and r Q = 0 
and is given by the first moment 

r=\ o rp(r)dr. (3.10) 

Table 3.1 shows the integration results for several values of /. 



Fraction of Largest 
Amplitudes 
Considered 


Mean Values 

(F y /VF 5 ") 


0.01 


2.359 


0.1 


1.800 


0.333 


1.416 


0.5 


1.256 


1.0 


0.886 



Table 3.1 Wave Amplitude Means 



c) Maximum Expected Wave Amplitude 

The expectation of the largest amplitude in a sample of N waves is found 
from the first moment of the probability distribution of the maximum amplitudes, r^. 
Results obtained from this computation are summarized in Table 3.2. 
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3. Empirical Formulation of Sea Spectra 

Several empirical formulas, based on the analysis of many wave records have 
been proposed to express the spectral density of the energy spectrum as a function of the 
wave frequency. The two most general cases use either wind speed, or significant wave 
height and significant wave period in the empirical formulas. The significant wave 
height and significant wave period, for most ocean engineering applications, is defined as 
the mean of the one-third highest waves and the mean of the wave periods associated 
with the one-third highest waves, respectively. 



Number of 
Waves 

N 


Maximum Wave 
Amplitudes 

('’max ) 


50 


2.12 


100 


2.28 


500 


2.61 


1,000 


2.78 


10,000 


3.13 


100,000 


3.47 



Table 3.2 Expected Maximum Amplitudes [Longuet-Higgins 1953] 



In the first case, the spectral density S( (d) is of the form 

CO 

where A and B are empirical constants, and V is the speed of the wind. 
In the second case S( (0) is given by 

r > 5 



(3.5) 



(3.6) 



where A and B are again empirical constants, T s is the significant period and H s the 
significant wave height. Table 3.3 presents a useful compilation of wave heights and 
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wave periods as a function of sea states and wind speed that may be used in the following 
formulas. 



Sea 


Wind 


Wave Height 


Wave Period 


State 


Velocity 


(feet) 


(sec) 




(kts) 










Average Significant H|/m 


T Tnax T au . 



0 


0 


0 


0 


0 


- 




- 




2 


0.05 


0.08 


0.10 


< 1.2 


0.7 


0.5 


1 


5 


0.18 


0.29 


0.37 


0. 4-2.8 


2.0 


1.4 




8.5 


0.6 


1.0 


1.2 


0.8-5. 0 


3.4 


2.4 


2 


10 


0.88 


1.4 


1.8 


1. 0-6.0 


4 


2.9 




12 


1.4 


2.2 


2.8 


l. 0-7.0 


4.8 


3.4 




13.5 


1.8 


2.9 


3.7 


1 4-7.6 


5.4 


3.9 


3 


14 


2.0 


3.3 


4.2 


1.5-7. 8 


5.6 


4.0 




16 


2.9 


4.6 


5.8 


2.0-8. 8 


6.5 


4.6 


4 


18 


3.8 


6.1 


7.8 


2.5-10.0 


7.2 


5.1 




19 


4.3 


6.9 


8.7 


2.8-10.6 


7.7 


5.4 


5 


20 


5.0 


8.0 


10 


3.0-11.1 


8.1 


5.7 




22 


6.4 


10 


13 


3.4-12.2 


8.9 


6.3 




24 


7.9 


12 


16 


3.7-13.5 


9.7 


6.8 


6 


24.5 


8.2 


13 


17 


3.8-13.6 


9.9 


7.0 




26 


9.6 


15 


20 


4.0-14.5 


10.5 


7.4 




28 


11 


18 


23 


4,5-15.5 


11.3 


7.9 



Table 3.3 Sea State and Wave Parameter Comparison [Berteaux 1976] 
a) Pierson-Moskowitz (PM) Formula 

A frequently used one parameter description of S( co) for a frilly developed 
sea, resulting from extended wind with unlimited fetch, is the PM spectrum. Using 
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Equation 3.5 as a basis, the single parameter used to describe the spectrum is the wind 
velocity V, with the empirical constants A and B having the values 0.008 lg 2 and 0.1 Ag 4 , 
respectively. The units associated with the gravitational constant g, and the wind velocity 
V must be the same for dimensional consistency. 



b) International Towing Tank Conference (ITTC) Formula 

Using statistical properties the relationship between the one-third 
significant wave height, Hm, and the wind speed, V, for the PM spectrum is 

V 2 

H m = 0.2092 — , (3.7) 

8 

therefore, the formula for the PM spectrum with the significant wave height as the single 
parameter is 

S(co) = 4e' (s//,,>4) . (3.8) 

co 

In this form, the empirical constants A and B have the values 0.008 Ig 2 and 0.0324g 2 , 
and the resulting one-parameter expression that uses significant wave height is referred to 
as the ITTC formula. 

Since this one-parameter formula describes a fully developed sea resulting 
from conditions that are rarely encountered (extended wind with unlimited fetch), the PM 
spectrum should be viewed as an asymptotic form. To overcome the limitations 
associated with the one-parameter spectral family, a two-parameter family, given by 
Equation 3.6, can be used. 



c) 



Bretschneider Formula 



The Bretschneider spectrum is a general form, and represents experimental 
data very well. The formula for the Bretschneider spectrum is 



S(o) - 42QQ ^ jl - 1050 /( 7 > 4 ) 
e 



(3-9) 



where H s and T s are the significant wave height and significant wave period, respectively. 
This expression can be used to represent developing, fully developed and decaying seas 
depending on the value of T s chosen, [Lewis 1989]. 
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d) International Ship Structure Congress (ISSC) Formula 

The ISSC spectrum is similar in form to the Bretschneider spectrum, with 
the exception that it is based on a mean frequency corresponding to the spectrum's center 
of area. This spectrum is represented by 



T a > 5 



(3.10) 



where H s and T are the significant wave height and mean wave period, respectively. A 
comparison of the ITTC, ISSC and Bretschneider spectrum for a sea state 3 condition is 
shown in Figure 3.2. 
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Figure 3.2 Comparison of Empirical Spectra 



e) Other Spectral Representations 

There are many other empirical relationships that have been used to model 
spectral characteristics of the sea, where each model is useful for a given area or sea 
condition that has specific characteristics. Examples of these are the Joint North Sea 
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Wave Project (JONSWAP) spectrum [Hasselmann 1973], the Ochi Six-Parameter Wave 
spectrum [Ochi 1076], the Wallops spectrum [Huang 1981] and the Generalized spectrum 
[Liu 1983]. The Wallops spectrum and the Generalized spectrum with variable 
exponents are adaptable for both deep and shallow water applications. 

It is common to use the recommended sea spectra from the ITTC and ISSC in ocean 
engineering design. For open sea conditions, the PM spectrum is recommended, and for 
fetch limited conditions either the Bretschneider, Ochi or JONSWAP spectrum is 
available. 

C. LINEAR REPRESENTATION OF SEA WAVES 

1. Spectral Approximations 

To properly estimate and cancel the wave induced disturbances acting on an 
underwater vehicle, some type of disturbance model is necessary. The necessity arises 
from the need to either embed the disturbance model in the vehicle's controller, or to use 
the disturbance model in a state estimator. It is known that a linear, Gauss-Markov 
stochastic process may be generated by sending white noise through an appropriate 
transfer function, where white noise is to mean a random process whose power spectral 
density is constant over the whole spectrum. Therefore, a linear approximation to the 
spectral representations in the previous sections can be obtained by sending a random 
white-noise signal through a second order filter, [Spanos 1981], This process can be 
written as 



In Equation 3.11, y(s) is the wave amplitude, q(s) is a white noise source with power 
spectrum 



yO) = h(s)q(s ) , 

with the linear approximation to the desired spectrum represented as 



(3.11) 



(D w (ry)=|/z(yty)| 2 cD w =||/iO-ty)|| 2 . 



(3.12) 




(3.13) 



and h(s) is a second order transfer function of the form 
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(3.14) 



h(s) = - 



2k(s/co o ) 



1 + 2 C(s/a) o )+(s/0) o ) 2 
Using equations 3.11-3.14 and substituting 5 = joy the approximate power 
spectrum can be given by 



<t> yy (co)=\h(jco)\ 2 <£> q(! (co) = - 



4k 2 (o)/ a ) ) : 



(3.15) 



(1 -(o)/o 0 ) 2 ) 2 +4 C 2 (co/o) 0 ) 2 
The power spectrum of this filter, <t>yy(Ci)), has zero energy for zero frequency and the 
maximum value of O xv (tyj) occurs at co ~ oy, for small values of £ This is desirable 
considering the shapes of the spectra displayed in Figure 3.2. The parameters C, and k in 
Equation 3.15, are found by minimizing the performance index 

J =j~($> yy Ua))-S(6))) 2 do). (3.16) 

2. State Space Formulation 



A linear state space model can be derived from equations 3.11 and 3.14. By 
defining the states as 



*1 =\[y(?)dr 
x 2 =y 



(3-17) 



the system states can be represented by the set of equations 



X, 




0 1 


"*1" 




0 


z= 


-co] -2 £co 0 _ 


+ 




*2. 




_* 2 _ 




_2kco 0 _ 



5 ’ = [0 




(3.18) 



The transfer function in Equation 3.14, is useful for representing the sea surface 
elevation, but it cannot be used to generate the wave velocity. This is seen by taking the 
limit 



limsy(s), (3.19) 

5 — 

as s tends to infinity. The result of this calculation is a constant, 2k&b, that cannot 
represent the oscillatory wave velocity. [Saelid 1983] has shown that using the proper 
transfer function 
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y 2£o(sla> 0 )s 

q (1 + 2 C(s/Q) o ) + (s/C0 o ) 2 ) 2 



(3.20) 



can solve this problem and [Riedel 1997] has demonstrated through the use of linear 
prediction theory that an eighth-order transfer function of the form 




q 



2^o\sJjo 2 )s_ 

\ + 2C{sIcq o ) + {sIcq o ) 2 \ ’ 



(3.21) 



allows extremely accurate matching of the target spectrum enabling fluid velocity 
prediction as much as one period ahead. 

Based on Equation 3.20, a generalized transfer function may be written in the 

form 



_ ks 2 
q (a + bs + s 2 ) 2 



(3.22) 



which implies 

(s 4 +2 bs 3 +(2 a+b 2 )s 2 +2 abs + a 2 )u f (s) = kq(s) . (3.23) 

[Astrom 1989] has shown that the parameters a, b and k, representing a PM spectrum 
based on sea state 3 conditions, can be found using Equation 3.16 to be a = k = 1 and 
b = 2. In this form. Equation 3.23 becomes an all-pole filter thus avoiding the problem of 
base period repetition of wave records, [Riedel 1997]. 

3. Spectral Modifications for a Moving Vehicle 



The spectrum that the vehicle “sees” while moving through a wave field with 
some forward velocity is not the same as that which a still vehicle would encounter. The 
actual spectrum that the vehicle encounters is a function of the vehicle’s forward speed U 
and its heading angle relative to the propagation direction of the sea waves, /?. The 
definition of the encounter angle (3 is given by 

p = n-(y-y/), (3.24) 

where y is the direction from which the waves are propagating referenced to the inertial 
reference frame, and y/is, the heading angle of the vehicle [Lewis 1988b], see Figure 3.3. 
The frequency modification that must occur to the disturbance spectrum is represented by 
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the Doppler equation, (Equation 2.56), given in the previous chapter. The Doppler 
equation, in an alternate more useful form is 

2 

co e ~co-— Ucosfi. (3.25) 

g 

Depending on the encounter angle, specific terms have been given to the 
orientation of the sea with reference to the vehicle, this is shown in Figure 3.4. As an 
example, consider the case where the encounter angle is zero, namely the vehicle is 
moving in the direction of wave propagation, this is referred to as a “following sea’\ In 
this case, as equation 3.25 indicates, the frequency of encounter, o\, is less than the actual 
wave frequency. For this case, it is interesting to note that the encounter frequency can 
be negative for large values of forward velocity U. 




Figure 3.3 Incident Wave Directions 

The methods for spectral approximation and state space realization of sea waves 
presented in this section are useful in simulation studies when the target spectrum is 
known, including the frequency shift. However, if the target spectrum is not known, 
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which is typically the case with a deployed vehicle, alternate techniques must be used. 
The following section presents one such technique which will allow a vehicle to estimate 
the encounter spectrum on-line. 



/ 3 = 90 ° 
Beam sea 




Figure 3.4 Encounter Angle 

D. WAVE MODELING USING RECURSIVE METHODS 

A typical experiment in system identification consists in recording a set of 
input/output data and fitting a parametric model. In discrete time, the attempt is to fit a 
Linear Difference Equation (LDE) model of the form 

y(t) + a { y(t - 1) + ... + a n y(t -n) = b x u(t -1) + ... + b n u(t -n) + £(0 , (3.26) 

with t : [0, °°] denoting the integer discrete time index, and e(t) an error term which 

accounts for the fact that the data never matches the model exactly. The problem is to 
estimate the parameter vector 

e = [a i ,...,a ll ,b l ,...,bj (3.27) 

from a set of data. This section will present the general procedure to estimate the 
parameters associated with a discrete linear model. 

Equation 3.26 may be written in regression form as 

y(t) = <p(t-l) r 0 + £(t) (3.28) 
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with 



<t>it-\) = [~ -y(t -n),u(t (3.29) 

being a sliding window of input and output data. Now the problem is to determine a 
technique to compute an estimate of the parameter vector 6 on the basis of the data set 
Z = {h(0),...,m(A 0, j(0),...,;y(AO}, if it is assumed that N data points are collected. 

If the problem is cast in a probabilistic framework, a probability density for the 
data set Z given the model 6 can be written as 



P(Z|0) = Pr(£(O), ...,£( AO) (3.30) 

with £{t) = y{t)-(j) T (t — \)6 . If we assume the disturbance term sequence to be 

Gaussian and white, then the density on the right hand side becomes 

Pr(£(0),...,£(A0) = f[Pr(£(0) (3-31) 

1=0 



with 



Pr(£) = 



1 

e 

27CCT 




As a consequence it can be seen that 



(3.32) 



. - i i r £|y(l)-/(,-l)*| 2 

Pr(Z|0) = Ce . (3.33) 



Minimizing the summation in the exponential can maximize this probability. So 
if the estimate of the parameters is defined as the vector that minimizes the probability, as 

6 = argmin e Pr(Z|0) (3.34) 



then it can be computed as a least squares solution 

M , 

§ = argmin e ]T|;y(O-0 r (O0|‘ • 



(3.35) 



The solution can be obtained using standard techniques, by writing 0 as the least 
squares solution of the system of equations 
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e 



(3.36) 



' y( 1) ' 




" m T ' 


y( 2) 


= 


m T 


_y(M)_ 




1 

/ — \ 

• j 

i 



or, compactly 



y = 0(9. 

The solution given by the pseudoinverse becomes 



d = (® T ®y i ® T y 



which can be written as 



(3.37) 

(3.38) 



1 M 

0 = — 2>(0*(0 r 



n-i 



1 M 



i = 1 



(3.39) 



This is true provided the error sequence e(t) is white and Gaussian. If the error is not 
white, then the estimate of the parameter vector d is biased. 

Now that we have set the foundation, the problem that arises is how do we 
estimate the parameters associated with a transfer function that will properly represent a 

A 

model of the seaway in a recursive fashion. If 0(t) is the estimate of the parameters at 

A A 

time t, the goal is to compute the successive estimate 6(t + 1) by updating 6(t) using the 
latest observations. It turns out that this problem can be put in a very nice framework that 
makes use of the considerations on the Kalman Filter. More on this approach will be 
discussed in Chapter IV. 

The auto-regressive (AR) model, with numerator equal to one, can be written as 

y(t) = W-l) T 0 + e(t) (3.40) 



with eft) a white noise sequence. If it is assumed that <?is constant. Equation 3.40 can be 
written in state space form, where the state is the parameter vector itself, as 

6 ft + l) = 0(t) 

(3.41) 
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This form is just a particular case of the stochastic state space model with A being the 
identity matrix, B=0 and C =0(t) T . Since this is a stochastic process, the Kalman Filter 



approach for the estimation of 0(t) is used, which leads to 



0(f + l) = 0(O + 



P{t)<Kt-\) 



X 2 +<p(t-l) T P(t)<p(t-l) 



(. y(t)-0(t) T ftt-i » 



F(t +1) =F( 



(3.42) 



X 2 +<p(t-l) T P(t)<p(t-1) 

with X 2 = E{|e(r)| 2 } Clearly, the parameter X is never known, and the need to know it 



can be eliminated by proper normalization. Dividing the numerator and denominator of 
Equation 3.42 by X 2 the following algorithm is obtained; 



3(t+i) = 3(t)+ 

p(t+i) = p(t)~ 



p{tm- 1) 



i+<p(t-i) T p(t)<p(t-i) 

p{tm - pm 

\+W-iy P(ty>{t-i) 






(3.43) 



which can be easily verified by setting P(t) = P(t)/ X } . This algorithm is called the 
Recursive Least Squares (RLS) algorithm, [Ljung 1987]. 

Using the above developed RLS algorithm the parameters of an AR model of the 

form 



w(t) + d } w(t -l) + ... + d N w(t - N) = e{t) , (3.44) 

with w(t) representing the sea surface elevation due to wave action, and e(t), a zero mean, 
white noise otherwise know as the innovation, may be developed. The main advantage of 
the AR model is the fact that the parameter estimation is a linear operation. In fact, if the 
numerator pertinent to the noise term is not one, then the error state model, Equation 
3.41, is not white and we cannot apply Kalman Filtering techniques directly. In the 
general case, the problem is nonlinear in the parameters and the Extended Kalman Filter 
must be applied. 

The problem that now arises is to determine the order (AO of the model to properly 
reflect the actual frequency spectrum of the time series w(t) while keeping the complexity 
of the model at a minimum. By computing the covariance of the innovation as the order 
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of the model is increased it can be shown that there is a value to which the covariance 
converges. Using surface wave elevation data obtained in Monterey Bay, from a 
Waverider measurement buoy, an Auto Regressive model of various orders was 
determined and the covariance of the innovation compared. The results of this analysis 
are shown is Figure 3.5. 

As can be seen, the covariance begins to flatten out around an eighth order model. 
This would have one believe that the correct order model to choose would be an eighth 
order model. However, as shown in Figure 3.6, an eighth order model fails to accurately 
reflect the actual spectrum As the order of the model is increased, the error in the 
matching of the actual spectrum decreases, however, not until a 100 th order AR model is 
computed, is the desired spectrum actually realized, see Figure 3.7. The issues associated 
with using a model order this high include parameter corruption due to noise causing 
inaccurate identification. 
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Figure 3.5 Comparison Of Innovation Covariance To AR Model Order 

There are those that will argue that the energy associated with the low and high 
frequency modes is minimal, and that proper modeling of those modes is not important. 
However, as was discussed in Chapter II, it is the low frequency modes of the wave train 
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that are considered shallow water, non-dispersive waves. It is these modes that will have 
the most effect on the submerged vehicle trying to maintain station. 
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Figure 3.6 8 th Order AR Model - Fit to Monterey Bay Data 
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Figure 3.7 100 th Order AR Model- Fit to Monterey Bay Data 
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Therefore, as in any design, trade-offs must be made based on analysis, as to how 
accurate the wave disturbance model must be compared to the resulting controller 
complexity. 

Using a model order much greater than that of an eighth order approximation 
creates difficulties in real-time embedded processes currently used on small underwater 
vehicles. However, using a stable, reduced order disturbance model in an estimator 
where measurements are available to correct the model errors provides a very good 
means of tracking and compensating for the disturbance, [Riedel 1998a], Further 
information regarding the estimation and accuracy of these models will be discussed in 
Chapter V. 

E. APPLICATION TO DISTRIBUTED SIMULATIONS 

This section will present a method by which realistic wave data as well as 
empirical relationships may be implemented into a distributed simulation. The purpose 
of this approach is to further the development of simulation capabilities for control 
system design, multiple vehicle coordination as well as mission planning. 

The simplest method of obtaining information about wave disturbances in a 
particular operating area is from a wave buoy. However, to use this information to 
simulate the disturbance forces and moments acting on a submerged vehicle, it is 
necessary to transform the wave elevation record into a water particle velocity record at 
the vehicle operating depth. There are two approaches to this problem. The first uses 
spectral analysis while the second uses Fourier analysis. 

Using the spectral analysis approach, the procedure is to first compute the power 
spectral density, <& nn ( (d), of the surface elevation, then modify the PSD by the appropriate 

depth related transfer function, |//(ry,Z)|“, and finally convert the new spectral density 

back to the time domain for replications of the subsurface water particle velocities. 
Although the resulting time series reflects the proper magnitude of the water particle 
velocities, the disadvantage to this method is that the phase information is lost, and the 
resulting subsurface velocity time series do not accurately reflect the motion caused by 
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the surface elevation times series. It is critical to ensure that the phase relationship 
between the horizontal and vertical wave induced velocities and accelerations are correct 
in order to provide realistic vehicle motion in simulation since, as Chapter II outlined, the 
6DOF EOM are coupled 

Through the use of the Fourier analysis method, the shortcomings of the spectral 
analysis are overcome and the disturbance phase relationship is maintained. In this 
method, a FFT of the surface wave record is taken, then the Fourier coefficients are 
multiplied by the appropriate value of the modifier for each frequency component. Once 
this has been completed, an inverse FFT is performed. The resulting disturbance records 
now reflect both the proper amplitude and phase relations. Figures 3.8-3.11 depict this 
Fourier analysis translation procedure for a wave elevation time series recorded in 
Monterey Bay, CA. 
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Figure 3.8 Wave Elevation Time Series, Monterey Bay April 1998 
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Figure 3.9 Surface Elevation PSD, 
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Figure 3.10 Sub-surface Water Particle Velocity Series, (H=45 m, Z=22.5 m) 
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Figure 3 1 1 Horizontal Water Particle Velocity PSD, O uu (/), (H=45 m, Z=22.5 m) 

The elevation data, displayed in Figure 3.8, was recorded by a Waverider® buoy 
in 45 meters (150 feet) of water. The resulting velocity record, Figure 3.10, is from a 
transformation of the elevation record for an operating depth of 22.8 meters (75 feet). 
Referring to the PSD of the transformed velocity record. Figure 3.11, it is interesting to 
note that in addition to the dominate peak frequency at 0 12 Hz , there are two additional 
lower frequency components at 0.05 Hz and 0.07 Hz. Recalling comments made with 
respect to shallow water waves in Chapter II, for this water depth only the frequency 
component less than 0.05 Hz may be treated as shallow water waves. 

To accurately reflect, in simulation, the wave field approximated by any of the 
empirical relationships discussed in Section B, Equation 3.3 as well as the relationships 
given in Table 2.2 are used. The stochastic nature of the sea waves is introduced in the 
simulation by including a random phase angle in the oscillatory 6 term given in Table 
2.2. The proper velocity amplitudes for this wave field are determined by using the 
vehicle’s depth in the exponential modifier. The correct phasing, with respect to the 
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vehicle, is obtained by using the vehicle's global position, X, in the 0 term associated 
with the water particle velocities. As a note, this approach models first order wave 
processes, however, to adequately represent the nonlinear second order wave effects, the 
Aruin Equation 3.3 is often varied to provide frequency bands of equal energy. 

F. SUMMARY 

This chapter has presented a detailed discussion of wave disturbances. It has 
described the statistical nature of the sea and presented several empirical relationships 
that may be used for approximating the spectral properties of the ocean. Methodologies 
for identifying and employing state space models of the sea, in control applications, have 
been highlighted. Approaches for using wave information, from empirical relationships 
or real data, in distributed simulations has been discussed. 
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IV. PARAMETER ESTIMATION 



A. INTRODUCTION 

As shown in Chapter II, accurate modeling of underwater vehicles has lead to the 
development of complicated equations of motion. Apart from the non-linearities inherent 
with underwater vehicle motion, the forces and moments acting on the submerged body 
are typically determined by a combination of theoretical and experimental results. 
[Healey 1992 and 1993] has suggested the use of three, independent, decoupled, 
equations of motion sets, which model speed, steering and diving control. The use of 
these simpler models, which describe only particular vehicle dynamics, is useful in 
control law design. 

Recent interest in underwater vehicle maneuvering and control in shallow water 
has generated the need for a greater understanding of vehicle dynamics in this regime. 
Specifically, improved vehicle models foster the development of sophisticated control 
architectures, which produce the high degree of autonomy necessary to allow vehicles to 
maintain acceptable performance in an ocean environment. Critical to the solution, since 
dynamic positioning of an underwater implies a nonlinear response, is the use of an 
adequate input-output mapping of the vehicle dynamics. 

In this chapter, a method for identifying the decoupled longitudinal surge motion 
dynamic parameters is presented. The identification is based on post-process Kalman 
filtering of data obtained from in-water vehicle experiments. Identification of the 
parameters associated with the surge equations of motion is performed, and a comparison 
between experimental data from in water measurements with the Phoenix AUV, and 
simulated results is conducted. Lastly, the nonlinear function that relates the commanded 
propeller speed to the required propulsion motor voltage is determined. This function is 
critical in the implementation of a real-time surge controller that will allow a vehicle to 
maintain position while disturbed by waves. 
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B. ESTIMATION METHODOLOGY 



Parameter estimation has been called a “can of worms” by Astrom 
[Astrom 1983], in which he refers to the difficulty of making theoretically sound 
methodology work with real data. With this said, many different techniques have been 
employed in the area of system identification or parameter estimation, see [Lung 1987, 
Gelb 1968 and Astrom 1989] for examples, to attack this daunting task. In this work, a 
recursive Kalman filter approach was chosen since this technique is suitable for real-time 
implementation. The Kalman filter method is similar to weighted least squares 
algorithms [Ljung 1987], and allows for the incorporation of system and measurement 
errors. 

In the Kalman filter algorithm, it is assumed that the parameter model is based on 
a nominally constant parameters set, where the state vector is the parameter vector 0, and 
the dynamics are described in discrete time by 



The system noise, w, and the measurement noise, v, are considered to be zero-mean, 
white noise sequences with associated covariance matrices Q and R, respectively. The 
recursive Kalman filter estimation equations as given by [Gelb 1974] are 



where K is a time varying optimal gain that produces a least squares solution for the 
parameter estimate, and P is the parameter estimation error covariance. These 
expressions are equivalent to the formulation by [Ljung 1987], where the “forgetting 
factor’ ‘ A, is related to the noise covariance R. This recursive algorithm is expressed as 



0(k + l)= 0 (k)+rw 

z(k+l) = h(k)0(k)+v ' 



(4.1) 




(4.2) 
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6{t) = d(t-\)+L(t)e{t) 

£{t) = [z{t)-h{t)d(t-l)} 

L(/) = P(t)3>(0 = P(f - 1 )$(o[A( 0/ +$ r (o^(f-i)®wl' 1 > ( 4 - 3 > 

p, t) = J_ \ P(t _ 1) _ P(t-Ymy* T (t)P(t-i) ~ 

ACO L KOZ+f^O^-lWO. 



where the standard Recursive Least Squares (RLS) algorithm is obtained for the special 
case of A = 1. 



C. PARAMETER IDENTIFICATION MODEL 

The surge model developed in Chapter II is in continuous time, however, 
real-time control and estimation is done in discrete time, therefore the set of equations 
given by Equation 2.91 must be converted to a digital form. To produce a digital set of 
equations, a standard Euler discretization can be used. Assuming a sampling period T, 
the discrete system of equations, disregarding the kinematic relation, for the surge model 
becomes 

u r (k + l) = (au r (lc)\u r (Jfc)| + F(k))T~ l +u r (k ) 

F(k + 1) = (— F(k) +?-u r (*)|n(*)| + £-n(k)\n(k)\)T~ l + F(k) ' 

The model presented in Equation 4.4 contains four unknown parameters, a, (3, y 
and T, that must be determined by experimental means. To properly identify the dynamic 
parameters associated with this mode, the model must be expressed in terms of the 
measurement variables only. The measured input-output data channels that are available 
for this identification are the relative longitudinal velocity, Ur(t), and the propeller 
revolutions, n(t). Since the propeller thrust value F, cannot be measured, the two 
first-order equations must be combined into a second-order equation containing only 
these two measured variables. Defining the change of variables 
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(4.5) 




T 



and combining the two first-order equations in Equation 4.4, a second-order discrete 
model in u r can be formed. 



u r (k + 2) = 2u r (k + l)—u r (lc) 

+ Ta[u r ( k + l)|n r (k + 1)| -u r (k)\u r (fc)|] 

+ Tb[u r (k)-u r (k + 1)] + T 2 ba[u r ( k)\u r (fc)|] ' (4 ' 6) 

+ T 2 c[u r (k)\n(k)\]+ T 2 d[n(k)\n(k)\] 

Unfortunately, the model in Equation 4.6 is nonlinear in parameters making it 
difficult to apply standard system identification techniques. However, the nonlinear 
model can be transformed to be linear in its parameters and variables, thereby allowing 
the use of well established estimation tools, by defining an additional change of variables 



as, 



and 



C,=Ta h l =u r (k + l)|tq (k + 1)| - u r ( k)\u r (fc)| 

C 2 =Tb h l -u r ( k ) -u r (k + 1) 

C 3 =T 2 ba h 3 =u r (k)\u r (k)\ , (4.7) 

C 4 =T 2 c h 4 =u r (k)\n(k)\ 

C 5 = T 2 d h 5 =n(k)\n(k)\ 



Z = u r (k + 2) — 2u r (k + 1) + u r (k) 

= >h h, h s ] . (4.8) 

e=[c, C 2 C, C. Cj 

With these definitions, the model can now be written in matrix notation as 

Z = h6. (4.9) 



The parameter vector, 6 contains five coefficients where only four parameters are of 
interest. The additional coefficient, C3, is a result of the nonlinearity associated with the 
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original model represented by Equation 4.9. The parameter vector, 6, will be estimated 
in a least squares sense using Equation 4.9, and in this case, one extra degree of freedom 
is present which may cause a slight decrease in the accuracy of the other parameters. 

In theory, the additional coefficient defines an implicit relationship between the 
other parameters of the model; in reality however, this additional coefficient lumps 
together any modeling errors and therefore it is not constrained thereby allowing 
convergence of the parameters in a least squares sense. The model parameters a, ft, y and 
T, are therefore determined from the four coefficients C/, C 5 , C 4 and C 2 , respectively. 

D. SYSTEM IDENTIFCATION 

1. Input Signal Design 

Prior to estimating the parameters associated with the surge model, data must be 
obtained for use in the estimation filter. This data must contain measured control input 
and response output variables from a “persistent’ ‘ excitation. It is a well-documented 
theorem that to ensure a unique, unbiased, least squares estimate, the system or plant 
must be persistently excited, [Astrom 1989]. In addition, the system identification should 
take place using open loop control, if possible, so that controller dynamics do not effect 
the results, [Ljung 1987]. 

Since the purpose of this piece of work is to determine the parameters of the surge 
dynamics model that will be the basis for a surge controller, the input data must 
persistently excite the vehicle over the expected frequency range of the surge velocity 
disturbance. For shallow water applications, the period of the surge disturbance that an 
underwater vehicle may encounter can range from approximately 4 to 40 seconds. It is 
necessary therefore, that the propeller revolution input to the vehicle, for identification 
purposes, also contains this frequency content. By selecting a square wave of various 
periods the control input was designed that contained the desired frequency components. 
A portion of the time series used for the control input is shown in Figure 4.1, and the 
frequency content of this input signal is depicted in Figure 4.2. 
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Figure 4.1 Control Input Time Series 




Figure 4.2 Control Input Frequency Content 
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As seen in Figure 4.2, the frequency components of this input signal range from 
0.02 Hz to 0.35 Hz. The dominant frequency was designed to be around 0.05 Hz to 
0.08 Hz,or a 12.5 to 20 second period. The purpose of designing the input signal in this 
manner is to ensure that the parameters were weighted in the range of the dominant surge 
period. With the input signal properly designed the system identification data runs can 
proceed. Note that considering the input to be n(t) or n(t)\n(t)\ makes no difference in the 
spectral content of the input signal. 

2. Data Collection 

A series of four in-water experiments were conducted in the Monterey Harbor 
Basin on March 5, 1999. The Phoenix vehicle was placed under active control in both 
heading and depth through the use of control surfaces, and the propeller RPMs were 
commanded through the use of an input file representing the input signal shown in Figure 
4.1. During these runs Phoenix carried its standard sensor suite, which includes a 
SonTek® ADV and a RDI® Navigator DVL. These sensors were used to measure the 
vehicle's response to the control input. Information concerning each of these sensors may 
be found in Appendix B. A sample of the measured input and output data obtained 
during run three is depicted in Figure 4.3. 

The input-output response shown in Figure 4.3 indicates that the vehicle is less 
efficient when operating astern. This is evident by observing that the magnitude of the 
astern velocity is significantly less than that of the forward velocity, for the same 
propeller revolutions. This decrease in efficiency is an issue that must be addressed 
during the estimation process. 

An additional item observed during data analysis is that the voltage required by 
the propulsion motors to obtain the same revolutions is different. The reason for the 
difference is attributed to the starboard propulsion train having to overcome a greater 
amount of friction. This friction is caused by misalignment that resulted from structural 
damage that was incurred by Phoenix early in its operating life. The voltage difference is 
shown in Figure 4.4. To account for the differences in propeller revolutions, an average 
input was used during the system identification. The speed from the two shafts was 
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Figure 4.3 Example of Measured Data, Propeller Revolutions, n, And Ground 

Velocity, u g , From The RDI 




Figure 4.4 Digital Voltage to Propulsion Motors 
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averaged, and the vehicle was treated as having a single input The parameter /? is then 
determined for the assumed single shaft system, with the recognition to the fact that the 
actual parameter for each shaft on the vehicle will be fill. 

Since the output of the digital controller is a commanded voltage, this voltage 
mismatch can be accounted for in the controller software. To ensure that the propellers 
produce approximately equal thrust, the function that maps the commanded revolutions 
from the controller to delivered voltage to the motors must be determined. Once this 
function is found, it can be coded into the vehicle control computer thereby ensuring the 
propeller mismatch is minimized. Another way of ensuring that the shaft speeds are the 
same is to “close the loop’ 4 around propeller speed. However, this solution adds 
additional system dynamics making the controller formulation more difficult. The 
methods used to determine this mapping will be discussed later in this chapter in 
Section E. 



3. Parameter Identification Results 

In the application of the parameter estimator, the measurement noise covariance, 
v, was set to 0.01, a constant scalar. The values of the diagonal elements, qa, of the 
model noise covariance matrix, Q, were chosen to match the bandwidth of the input 
signal With this choice, a tradeoff between convergence, stability and precision of the 
estimates is made. The following results are presented for one of the many data runs and 
selected values of qu that produced the “best” results in parameter estimation. The term 
“best” is a subjective measure based on a balance between the whiteness of the residuals, 
the values of the parameters compared to parameters previously identified and a 
comparison between measured data and simulated results. The final results of the 
selected parameters from the estimation process are shown in Table 4.1, including 
statistics on the estimation error and the residual e. 

The evolution of the parameter estimates and the diagonal elements of the 
covariance matrix, P, associated with these parameters is shown in Figures 4.5 and 4.6, 
respectively. Referring to Figure 4.5, it can be seen that the parameters do converge, but 
exhibit some slight fluctuations. The noise observed in the evolution of the covariance 
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Table 4.1 Parameter Estimation Values and Residual Statistics 
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Figure 4.5 Parameter Evolution 
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Figure 4.6 Covariance Evolution 

matrix, Figure 4.6, is caused by the rather large values of the diagonal elements of Q, 
qu = 5.0. The large values needed for these were necessitated by the additional degree of 
freedom in the filter, the uncertainty associated with the thrust reduction term and the 
large bandwidth of the input signal. 

The performance of the filter can be determined by computing the autocorrelation 
of the residuals. The autocorrelation provides a measure of how the value of a random 
variable at a time t, will influence its value at some future time, t+T. If the filter is 
properly tuned, the residuals should be white. For the residuals to be classified as white, 
there will be ideally zero correlation at any time shift T, and the time series will be highly 
correlated at r= 0. Figure 4.7 indicates that the residuals are not white, but exhibit some 
correlation. The non-whiteness of the residuals indicates the lack of modeling capability. 
Since Kalman filter theory assumes that any measurement noise has a Gaussian 
distribution, the measurement noise model is corrupted by unmodeled, colored noise. 
This corruption of the measurement model caused the residuals also to display colored 
noise properties. Also, since the reduction in propeller efficiency when the vehicle 
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operates astern is not known, additional modeling errors are introduced which effects 
filter performance. 

Despite the issues associated with the non-whiteness of the filter residuals, 
apparent satisfactory parameter identification was obtained. The parameter estimation 
was deemed satisfactory by comparing the measured vehicle response to simulated 
response results. Using the identified parameters, listed in Table 4.1, and the measured 
propeller input shown, in Figure 4.3, a simulation was conducted with the results 
presented in Figure 4.8. 

The response traces shown in Figure 4.8 indicate that the identified parameter set 
provides a reasonable predictive response. Recalling that the parameters are estimated in 
a least squares sense, the errors between the two traces are acceptable. The two 
responses are not lagged indicating good agreement in the propulsion system time 
constant and vehicle drag coefficient. The error in response amplitude is attributed to the 
uncertainty associated with the identification of the /? and /parameters. 




0 500 1000 1500 2000 2500 3000 

Figure 4.7 Autocorrelation of Estimation Filter Residuals 
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Figure 4.8 Comparison of Measured and Simulated Vehicle Response. 

The uncertainty in each parameter may be estimated by using the magnitude of 
the covariance matrix diagonal corresponding to the estimated parameter. The calculated 
uncertainty for each of the parameters, based on the covariance levels, is shown in 
Table 4.1. Although the level of uncertainty with these parameters may seem excessive, 
it is not too different from the “standard” level of uncertainty associated with underwater 
vehicle parameter estimates, which typically is around ±40 %. It is this variability in 
parameters, as well as other items, that produces the need for robust control law design, 
which will be discussed in the following chapter. 

E. PROPULSION SYSTEM BALANCING 

As outlined earlier in this chapter, the required propulsion motor voltage to 
produce a desired propeller revolution is different for each shaft. Under position control, 
this difference if not accounted for could cause unsatisfactory results since for position 
control the loop is closed around position error and not propulsion shaft speed. 
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In the control system architecture of the Phoenix AUV, the shaft speed is 
controlled by a motor voltage command which originates from the execution level 
computer. The particular control law that has been implemented in the vehicle calculates 
this voltage; therefore, the relationship that maps desired shaft revolutions to voltage 
must be determined. 

By fitting a curve to an overlay of data obtained during the four, system 
identification runs, the relationship between propeller “revs” and motor voltage for each 
shaft can be determined. Using the MATLAB® polyflt algorithm it was determined that 
a third order polynomial adequately mapped the two input-output relationships. The 
nonlinear functions for these relationships are, 

v„.,„ =0.022607n>„ + 0.000194^,, + 0.832982«„, 

V„ f = 0.022356,4, - 0.009507«1 + 0.345086n„_ 

with the graphical results of the “goodness of fit” shown in Figure 4.9. 





Shaft Revs (rps) 



Figure 4.9 Nonlinear Shaft Revolution to Motor Voltage Function 



As a note, the polynomial curves in Figure 4.9, indicate that the starboard shaft 
will produce approximately 8.5 rps for an applied voltage of 24 volts, while the port shaft 
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produces almost 10 rps for the same applied voltage. Although there is no data displayed 
in Figure 4.9 to validate this curve at high speed, the shaft revolution values estimated by 
these curves have been observed during full power trials. 

F. SUMMARY 

This chapter has presented the work performed to identify the parameters 
associated with the longitudinal dynamics of the NPS Phoenix AUV from open water 
experiments. The results show the necessity of a non-linear dynamics model and of the 
need to include a force lag and a thrust reduction term in the propulsion model. The 
results obtained are satisfactory although the filter residuals appear to be non-white due 
to limited modeling capabilities. Regardless of the level of whiteness, it must be stressed 
that the simple model identified is very useful for control law development. 

One interesting point observed during the estimation process was that taking Q as 
a diagonal matrix with equal values produced better results in the sense that the residuals 
were minimized and the filter produces stable parameter evolution. This reflected a 
priori, equal certainty in each initial parameter estimate. Moreover, the interrelationships 
between the parameters implied a uniform characterization of the modeling noise in the 
filter. The tuning of the estimation filter proved to be an important step towards 
obtaining a good model. 

Finally, the methodology presented here can be extended to other types of vehicle 
motions and represents a basis for the development of models and identification 
techniques that can be applied to other coupled or decoupled motion models. 
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V. DISTURBANCE REJECTION THEORY 



A. INTRODUCTION 

Many different techniques have been applied to the disturbance rejection problem 
in the past. For years, passive techniques employing improved electro-mechanical design 
were used. With current advancements in the computing industry, adaptive methods have 
become a popular means of active control. Recent research in the area of nonlinear 
control has shown that Variable Structure or Sliding Mode Control provides a very robust 
method of disturbance compensation. 

This chapter will begin with an overview of disturbance rejection theory 
developed to date. This will provide the reader a foundation from which to see the 
extensions made here in regards to the development of disturbance rejection techniques 
for small underwater vehicles. 

Next, a series of three case studies will be presented. These case studies will 
demonstrate the major approaches available in the design of disturbance rejection 
compensators, including the effect that the spectral content of the input disturbance has 
on controller performance. 

Lastly, the results of a simulation study, for the design of a station keeping 
controller that will be used on the NPS Phoenix AUV, will be discussed. 

B. OVERVIEW OF CLASSICAL CONTROL TECHNIQUES 

The most intuitive and oldest means of eliminating the effects of a disturbance is 
to attempt to attenuate the disturbance at the source. This often translates to corrective 
measures in the system For example, modifying the electronics in a sensor so that the 
noise is reduced, is one common application of this technique. Other examples are 
reducing friction forces in a servo by using better bearings, or moving a sensor to a 
position where the disturbances are smaller. Although this method of reduction at the 
source is beautiful in its simplicity, it is often impossible to achieve. 
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Feedback Control 



If the disturbances cannot be rejected at the source, feedback control can be used. 
For this method, the manner in which the disturbance enters the system must be known, 
and the system must be both controllable and observable. In this way, the effects of the 
disturbance can be mitigated by using local feedback. 

The classical control problem. Figure 5.1, simply stated, is, given a plant model 
G p , design a controller, G c> , such that the closed loop system 

1. is stable and exhibits some level of robustness against plant parameter variations; 

2. accurately tracks the reference input signal, r; and 

3. rejects the disturbance d, and noise v. 

The solution to this problem is accomplished, in general, by selecting a controller such 
that a high loop gain is obtained over the frequency range of the input disturbance, while 
obtaining a low loop gain over the range of the frequencies associated with the 
measurement noise. 




Figure 5. 1 Block Diagram of a Feedback Controller 

A proper design can be obtained by simultaneously solving the two performance 
criterion. 
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(5.1) 




over the range of frequencies of interest. This approach is typically referred to as "loop- 
shaping,” and any of the standard control design techniques can be used to accomplish 
this goal, [Ogata 1990]. 

To eliminate any step or bias error (or to follow a step input), an integrator is 
added to the system The integrator directly changes the system's sensitivity to 
disturbances, specifically targeting bias disturbances. Integral control, in application, 
needs to be managed carefully using anti-windup methods. This compensation technique 
does not specifically reject the disturbance, but acts to alter the resulting dynamics of the 
system when a disturbance is present. 

2. Feedforward of Directly Measured Disturbance 

If a disturbance can be measured or estimated, feedforward control is a useful 
method of canceling its effects on the system response. Unlike feedback control, this 
method is advantageous in that it is implemented by approximately compensating for 
disturbances before they are sensed. In feedforward control, a signal from a measurable 
disturbance maybe used to generate an appropriate control force to counteract or mitigate 
the effects of the disturbance. It minimizes the magnitude of the output for the 
disturbance input without the use of error integration. 

Feedforward control alone can minimize transient errors but there are no 
guarantees of its accuracy due to its open-loop nature. Thus, feedforward alone is 
unrealistic for most applications with unsuitable open-loop dynamics. For this reason, 
feedback control is often used together with feedforward control to compensate for the 
latter's inaccuracies in minimizing the error. This approach is shown in Figure 5.2. 
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Figure 5.2 Feedforward-Feedback Controller Block Diagram 



C. MODEL BASED CONTROL 

The state-space representation of a linear system may be defined by the following 
set of equations: 

x(t) = Ax(t) + Bu(t ) 

y(t) = Cx{t) + Du{tY ' ’ 

where A, B, C and D are determined from the differential equations describing a system. 
A block diagram depicting this system is shown in Figure 5.3. 
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1. Linear Quadratic Regulator Control 

Performance indexes are a way of obtaining desirable output regulation without 
requiring excessive input signals. One of the most common and useful is the LQR 
performance index defined by 

'/ 

J lqr = x T Mx + j{x T Qx+u T Ru]dr (5.3) 

o 

where M and Q are real, symmetric positive semi-definite matrices and R is a real, 
symmetric positive definite matrix. M is called the terminal penalty matrix, Q is the state 
weighting matrix, and R is the control weighting matrix. 

Using a performance index of this form, subject to the system dynamics given in 
Equation 5.2, a linear state feedback (LSF) control law of the form 

u(t) = Kx(t) (5.4) 

was found to be optimal for minimizing Jlqr- This control law results in the block 
diagram depicted in Figure 5.4. 




Figure 5.4 Block Diagram of Closed-Loop State-Space System with LSF 

If a unique positive definite solution to the steady-state matrix Riccati equation, 

Q-PBR- i B t P + PA + A t P = 0 , (5.5) 
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represented by P, existed, then the LSF control law in Equation 5.4 results in the 
minimization of Jlqr • The minimization of Jlqr implies a desire to minimize both 
excessive output excursions and the control effort required to prevent such excursions. 
The adjustable weights M, Q, and R can be used to obtain an appropriate compromise 
between these two conflicting goals. The optimal control for this problem then becomes 

u(t) = -R' I B T Px(t). (5.6) 

Use of LQR control results in the optimal gain K and optimal pole positions. This 
method works for time-invariant or time-varying systems and is just as easy for MIMO 
systems as for SISO systems. Like the classical control methods previously discussed 
however, this control method is a feedback approach which attempts to tailor the system 
response while not specifically rejecting a disturbance. However, if the state dynamics 
matrix is augmented, to include the internal states of the disturbance signal, the 
disturbance can be rejected. 

2. LQR Control with Disturbance Feedforward 

Given the state-space equation of the system, 

x = Ax+Bu + F d d (5.7) 

if B and F are collinear, then the system can be rewritten as 

x = Ax + B(u + ad), (5.8) 

and the control law can be written in the form 

u=-Kx-od, (5.9) 

as long as the disturbance is measurable. 

However, if B and F are not collinear, then direct feedforward cannot be used and 
a disturbance estimator must be employed. With this approach, the state space system 
must be augmented with a disturbance model, and a controller designed based on this 
augmented system. With a design of this form, the LQR controller can account for the 
effects of the unwanted input. This can be done only when some assumptions about the 
form of the disturbance model can be made. 

The disturbance state z, with internal dynamics Ad, may be represented by 

z = A d z, (5.10) 
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where 



d = C d z. (5.11) 

Augmenting the system states, x, with the disturbance states, z, yields the following state 
space equation 

which may be expressed compactly as 

x a (t) = A a x a (t) + B a u(t)+w, u = -K l x-K 2 z , (5.13) 

with Ki and Kz representing the feedback and feedforward gains, respectively. 

With the disturbance dynamics included in the state dynamics matrix, the cost 
function, Jlqr, is minimized to determine a new set of gains based on the augmented 
system. A disturbance state estimator is necessary to provide the internal disturbance 
states and the LQR controller uses these estimated states in the state feedback loop. 
Figure 5.5 depicts this system. Note, that the augmentation of the states has no effect on 
the dynamics of the disturbance estimator since it is uncontrollable from u(t). The 
disturbance state estimation is driven by measurements from the output or disturbance as 
available. With the formulation shown in Equation 5.13 the effects of the disturbance can 
be reduced, and in theory cancelled if a perfect disturbance model is available. 

To this point, the control methods discussed have been based on linear or 
linearized systems. The ability to use these tools for the design of controllers that will be 
used on nonlinear systems is some what limited. In general, the linear systems may not 
very robust to model mismatch which can result in system instabilities, although, 
robustness measures can be implemented into the design process by HVTL and LMI 
techniques, [Silvestre 1998a,b]. 

3. Nonlinear Methods 

In nonlinear control, the concept of feedback plays a fundamental role in 
controller design, as it does in linear control. However, the importance of feedforward is 
much more marked than in linear control. Feedforward is used to cancel the effects of 
known disturbances and provide anticipatory corrections in tracking behavior. Very 
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often, it is impossible to control a nonlinear system without feedforward compensation. 
Note that a model of the plant is always required for feedforward compensation, although 
this model need not be very accurate. 




Figure 5.5 Block Diagram of Closed-Loop State-Space System with LSF and Estimated 

Disturbance Feedforward 

There is no general method for the design of nonlinear controllers. What is 
available to the designer is a collection of tools that are applicable to particular classes of 
nonlinear control problems. These nonlinear design tools can be placed in one of five 
categories: Trial-and-error, Feedback/Input Linearization, Adaptive Control, Robust or 
Sliding Mode Control and Gain-scheduling. Unlike the linear control discussion, the 
sections that follow will briefly discuss some of the salient points of these techniques. 
For further detailed information on each of these design tools, the reader is referred to 
[Slotine and Li 1991]. 
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Trial-and-error can be used to develop controllers. This method is similar in 
approach to linear lead-lag compensator design using Bode plots. The goal is to use 
analysis tools, (i.e., phase plane, describing function, Lyapunov analysis), to assist in the 
search for a control solution that can be qualified by analysis and simulation. Experience 
plays a major role in this technique, which, for complex system typically fails. 

Feedback linearization deals with techniques for transforming complex models 
into equivalent models of a simpler form, [Slotine and Li 1991]. In this nonlinear design 
methodology, the idea is to first transform the nonlinear system into a full or partially 
linear system. Once this has been done, then any of the linear design tools may be used 
to develop the necessary control system. Two draw backs of this method are that it 
typically requires full state feedback and it is not very robust to parametric uncertainty or 
disturbances. These draw backs can be overcome by the use of either robust or adaptive 
control methods. 

For uncertain or time varying systems adaptive control is very useful. Current 
adaptive control design applies mainly to systems that have well known dynamics, but 
unknown or slowly varying parameters. Adaptive controllers, whether developed for 
linear or nonlinear systems are inherently nonlinear. For nonlinear systems, adaptive 
control can be viewed as an alternative to robust nonlinear control. 

Robust nonlinear control techniques have proven very effective in a variety of 
practical control problems, [Healey 1993, Marco 1996, Yoeger 1991, Young 1996]. The 
controller is designed based on consideration of both the nominal model, and some 
characterization of the uncertainties associated with the model. Sliding Mode Control 
provides a systematic approach to the problem of maintaining stability and performance 
in the presence of modeling inprecisions. 

Gain scheduling is an attempt to apply linear control methods to the control of 
nonlinear systems. It was originally developed by the aircraft industry for the control of 
high precision aircraft. The idea is to select a number of typical operating points which 
cover the systems range of operation. The plant is then linearized and a controller 
designed for each of these points. Between operating points, the gains of the 
compensator are scheduled resulting in a global controller. The main problems with gain 
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scheduling is that it has only limited stability guarantees for nonlinear operations, and the 
computational burden of computing many linear controllers. 

D. CASE STUDIES IN DISTURBANCE REJECTION 

This section will present and discuss three distinct cases of disturbance rejection 
for three different disturbance inputs. The performance of each of the control designs 
will be evaluated by using the nonlinear EOM, Equation 2.94, in simulation studies. 
During the simulations, each control design will be subjected to a simple harmonic 
disturbance input, a PM spectrum based disturbance input and real disturbance data 
obtained from Monterey Bay. Sample data records for each of these disturbance inputs 
are shown in Figure 5.6. 

The three cases are summarized below: 

• Case I: High Gain LQR Control. Equation 2.94 will be linearized around a 
nominal operating point. Based on this linearized model the control gains for 
a full state feedback controller are calculated using a LQR method. 

• Case II: LQR Control With Estimated Disturbance Feedforward. Employing 
the linearized model from Case I, the system is augmented with an AR model 
representing the disturbance dynamics. The augmented system is used to 
calculate the control gains, and the AR model is used as a basis for a 
disturbance estimator. This controller uses full state feedback, with estimated 
disturbance state feedforward in the control calculation. 

• Case HI: Sliding Mode Control (SMC) With Measured Disturbance 

Feedforward. A model based sliding mode controller will be developed which 
embeds the disturbance in the control formulation. This controller relies on 
full state feedback with measured disturbance feedforward. 

1. Case I. High Gain LQR Control 

Using the 1 DOF surge EOM as a model, it is necessary to linearize this system of 
equations in order to use linear techniques in the controller design. Linearization can be 
performed a with a variety of approaches, stochastic [Leira 1987], harmonic 
[Heyns 1995] or nominal operating, pointwise linearization, condition [Riedel 1998a], 
but for this design, since the control law will attempt to allow an AUV to hold position 
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Figure 5.6 Sample Disturbance Input Time Series 



with zero ground velocity, the system is linearized about the steady state solution to 
Equation 2.94 while in the presence of a steady current. Performing this steady state 
analysis yields, 
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as the nominal operating point, where n Q must be real and the same sign as u r , 0 . Using a 
standard Taylor series linearization, the linearized system of equations, in state space 
from, can then be written as 



X 




0 


1 


0 


X 




0 




T 


u r 


= 


0 


2asgn(u ro )u ro 


1 


U r 


+ 


0 

7 P 

-u ro s ign(n 0 ) + 2 — sgn(n 0 )n 0 
_T T 


n + 


0 


F 




0 


fw 


-1 
T . 


F 






0 





~i 


0 


o' 


X 


y = 


0 


1 


0 






0 


0 


1 


F 



(5.15) 

If it is assumed that the vehicle will be operating in a -0. 1 m/s steady current, and the 
parameters a, ft, /and rare available then Equation 5.15 can be evaluated numerically. 

Since the parameters identified in Chapter IV were obtained from a discrete filter, 
and it is desired to implement the to be developed control law in a digital computer, the 
state space equations must be converted into a discrete form. Using a standard Euler 
discretization, Equation 5.15 can be represented in discrete form as 
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Using standard optimal control techniques the solution for the optimal (LQR) 
controller can be found as 
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n = -R~'B T Sx 



(5.17) 



where S is found by solving the steady state algebraic Riccati equation (ARE) 

A t S + SA-SBR-'B t S+Q=0 , (5.18) 

for the positive definite matrix S. In Equation 5.26, Q is the weighting matrix on the state 
error, and R is the weighting scalar, since this is a single input system, that invokes a 
penalty against the control effort. The LQR approach will always yield a stable system, 
as long as the Riccati equation provides a positive definite solution matrix S, for which 
the system must be controllable and full state feedback available. 

The following sections will show the simulated performance of the LQR 
controller when subjected to various disturbance inputs. The purpose of the simulations 
in these sub-cases is to provide a baseline by which to compare the performance the 
controllers developed in Cases II and III. 

a) Monochromatic Disturbance Input ( Case la ) 

Using the sine wave disturbance input depicted in Figure 5.6, the LQR 
controller was formulated and simulated for various control weighting values (/?). A plot 
of the position response for one of these simulations is shown is Figure 5.7. This 
response, the "best" of the many simulations, is the result of a high gain controller. As 
can be seen, the oscillations about the commanded position (0 meters) are significant and 
poor disturbance rejection is obtained. In addition, there is an obvious offset caused by 
the mean disturbance. This offset can be corrected by incorporating integral control into 
the LQR design, however, integral control will not correct the severe positional 
oscillations. 

During these simulations, the propellers were not limited, i.e., no 
saturation. The control input necessary to obtain the positioning shown in Figure 5.7 is 
displayed in Figure 5.8. As shown is this figure, the propeller oscillations are extreme 
considering that the model and parameters used in the simulations are based on the NPS 
Phoenix AUV which has a maximum propeller revolution of 800 rpm. Once again, it 
must be pointed out that the purpose of the studies in Case I is to obtain a baseline for 
comparison. 
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Figure 5.8 Propeller Response for Case la with Monochromatic Disturbance Input 
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b ) PM Spectrum Based Disturbance Input 

Using the same three state model as with Case la, simulation studies using 
a PM spectrum based disturbance input were conducted. The input disturbance was 
based on a significant wave height of 1 meter in a water depth of 45 meters. The vehicle 
was assumed to be operating at a 25 meter depth. The goal of this simulation study was 
to determine the control performance based on a disturbance input which contained a 
range of frequencies which the vehicle may encounter. 

Using the controller design resulting in the responses displayed in Figures 
5.7 and 5.8, a simulation was conducted resulting in the position response shown in 
Figure 5.9. In this particular simulation, the standard deviation of the position response is 
significantly reduced due to the magnitude of the disturbance input. Comparing 
disturbance inputs between Cases la and lb, it can be seen that the magnitude of the PM 




Figure 5.9 Position Response for Case lb with PM Spectrum Based Disturbance Input 

based disturbance is 15 times less than that of the monochromatic disturbance. This 
reduction in oscillation magnitude is also reflected when comparing position responses. 
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The propeller input response shown in Figure 5.10, is also significantly less, but still 
exceeds the maximum propulsion system input of 800 rpm. 
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Figure 5.10 Propeller Response for Case lb with PM Based Disturbance Input 
c) Monterey Bay Disturbance Input 

Using the control design that resulted in the responses displayed in Figures 
5.7-5.10, a third set of simulations was conducted. In this set of simulations, transformed 
wave buoy data obtained in Monterey Bay, CA was used as the disturbance input. This 
data was obtained from a Datawell® Waverider Buoy deployed April 9, 1998, from the 
research vessel R/V POINT SUR, during an NPS oceanography class (OC4610) cruise, 
under the direction of Prof. Thomas Herbers. 

The wave buoy, according to Defense Mapping Agency navigation charts, 
was deployed in approximately 45 meters of water. The wave elevation data obtained 
from the buoy was transformed to a subsurface velocity record, at a depth of 25 meters, 
using the procedure outlined in Chapter IV. A sample of the resulting time series, with a 
-0.1 m/s steady current superimposed, was displayed in Figure 5.6 
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The position response using this input disturbance is shown in Figure 5.11. 
The standard deviation of this response is approximately one-half of the response 
obtained from the PM based disturbance input with the same control design. Once again, 
this is due to the fact that the standard deviation of the Monterey Bay input disturbance is 
about one-half the PM disturbance. 
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Figure 5.11 Position Response for Case Ic with Monterey Bay Disturbance Input 

It is interesting to note, that although the position response reduced by a 
factor of two, when compared to the PM based case, the propeller input response did not, 
see Figure 5.12. This is due to the fact that the frequency content of the input 
disturbances is much different. This is evident by referring back to Figure 5.6. The Bay 
data contains more high frequency components causing the propulsion system to respond 
much more. 

Since the propulsion system response is still in excess of maximum output, 
tuning of the controller gains must be performed to bring the propeller rpms within limits. 
By adjusting the input weighting scalar R, and reducing the controller gains the maximum 
commanded propeller revolutions can be reduced as well as reducing the sensitivity of 
the controller to high frequency "noise.” This reduction of propeller input is at the 
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expense of increased position error. These results are shown in Figures 5.13 and 5.14. 
As shown in Figure 5.13, the standard deviation has increased by a factor of two in order 
to keep the propeller revolutions within propulsion system limits. 
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Figure 5.12 Propeller Response for Case Ic, Monterey Bay Disturbance Input 
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Figure 5.13 Position Response for Case Ic, Monterey Bay Disturbance Input 
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Figure 5.14 Propeller Response for Case Ic, Monterey Bay Disturbance Input, rpms 

Within Design Limits 

Figure 5.15 shows graphically the relationship between the level of control 
input and the level of disturbance rejection for the standard LQR solution subjected to 
real wave data from Monterey Bay. The position covariance is normalized by the 
covariance of the "free floating" or uncontrolled response of the vehicle, and the input 
covariance is normalized by the maximum rpm available from the propellers. This 
analysis can give a "feel" for how tight a control law must be provided to achieve a 
reasonable disturbance rejection. 

2. Case II. LQR Control with Disturbance Estimation Feedforward 

The problem that now must be addressed is how to achieve better performance. It 
has been shown , that by embedding an estimator of the disturbance into the control 
system design, improved performance may be obtained [Grimble 1995, Riedel 1998a]. 

As outlined in Chapter III, an AR model of the wave disturbance may be written 
in state space form as 

X w (k + l) = A w X w (k) + B w v(k) 

y w (k) = u f (k) = C w X w (k ) ’ ( ' } 



107 




Figure 5. 15 Comparison Of Control Input Covariance To Normalized Vehicle Position 

Covariance, Monterey Bay Wave Data 

Augmenting the vehicle state equations with the disturbance state equations, a 
new control law may be developed using the estimated disturbance states. Defining the 
new state vector as 

X** = 

where the disturbance states are given as 

X w = K(* + Af-l),...,X w (* + l)r , (5.21) 

the new control law may be designed, using the separation principle, assuming all states 
are measurable. As in the previous optimal control discussion (Case I), the ARE is 
solved to obtain the appropriate gains for the augmented system. This augmented system 
is represented by 

v(k). (5.22) 



X aug (k + l) = 



x vehicle 

0 



FC« 

A, 



*--(*) + 



‘ aug 



B 



vehicle 

0 



n(k) + 



v vehicle 

X... 



(5.20) 
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With the control law determined, the estimator must be designed. Using optimal 
estimation theory, an estimator of the form 

X w (k + 1) = (A w - LC W )X w (k) + Lu f (k) , (5.23) 

where u/k) is the current disturbance measurement, is developed. This estimator is used 
in conjunction with the control law developed, and its implementation, in block diagram 
form, is represented by Figure 5.5. 

a) Monochromatic Disturbance Input 

To display how well this design procedure can work if an accurate model 
of the disturbance is available, consider the case of the monochromatic input disturbance. 
Since the precise model of this disturbance is known, when this model is embedded in the 
control system design, perfect cancellation of the wave disturbance effects on vehicle 
positioning may be obtained. These results are displayed in Figures 5.16 and 5.17. Now 
these results are for demonstration purposes only, and perfect cancellation of the wave 
disturbances will not be possible since an exact model of the random sea is not available. 
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Figure 5.16 Position Response for Case Ha, Monochromatic Disturbance Input 
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Figure 5.17 Propeller Response Case Ha, Monochromatic Disturbance Input 

Although the propeller response is almost identical to the results displayed 
in Case la (Figure 5.8), by having an estimate of the disturbance states to feedforward the 
propeller input is properly phased to cancel the disturbance. 

b) PM Spectrum Based Disturbance Input 

Using this approach of a model based disturbance estimator with a LQR 
controller appears to be an excellent method of canceling the disturbances acting on an 
underwater vehicle, that is if the model of the disturbance is known. If the exact model 
of the disturbance is not known, the question is; Is improved disturbance rejection with 
this method possible? 

Adopting the AR modeling techniques presented in Chapter IV, a sixth 
order AR model for the PM based disturbance was developed. Using this linear model of 
the disturbance dynamics and the same input weighting scalar as was used in Case lb, a 
combined controller/estimator was developed. Using this developed compensator in 
simulation, improved performance was observed, see Figure 5.18. The position response 
with the augmented disturbance model improves by a factor of 1.5. 
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Figure 5.18 Position Response for Case lib, PM Based Disturbance Input 

The improvement in propeller input follows the trend displayed in the 
comparison between Cases la and Ha. The standard deviation of the input response has 
not changed significantly, as evident in comparisons between Figures 5.19 and 5.10. 
What has changed is the control input phasing, again due to the disturbance feedforward, 
thus allowing this design method, even with a low order disturbance model, to obtain 
improved disturbance rejection. 

c) Monterey Bay Disturbance Input 

Using identical weighting values that went into the design of the control 
laws used in the simulations presented in Case Ic, Figures 5.11-5.14, and a sixth-order 
AR model representing the Monterey Bay disturbance, improved performance was again 
realized. As can be seen in Figure 5.20, there is a 150% improvement in station keeping 
as compared Case Ic, and the control input requirements are significantly less, see Figure 
5.21. Although the standard deviation of the commanded control input is well within the 
maximum revolutions able to be provided by the propulsion system, there are some 
inputs which exceed the limit of 800 rpm. In order to bring the commanded control input 
within limits, as before, the control gain must be reduced. 
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Figure 5.19 Propeller Response for Case lib, PM Based Disturbance Input 
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Figure 5.20 Position Response for Case lie, Monterey Bay Disturbance Input 

With the control gains adjusted so that the commanded control input remained 
within propulsion system limits, the positional error increased by a factor of two, while 
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the control input reduced by a factor of three The results of this tuning are shown in 
Figures 5.22 and 5.23. This reduction in control effort is particularly important given the 
fact that power consumption is the downfall of AUVs. 
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Figure 5.21 Propeller Response for Case lie, Monterey Bay Disturbance Input 
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Figure 5.22 Position Response for Case lie, Monterey Bay Disturbance Input, rpms 

Within Design Limits 
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Figure 5.23 Propeller Response for Case lie, Monterey Bay Disturbance Input, rpms 

Within Design Limits 

3. Case HI. Sliding Mode Control with Measured Disturbance 
Feedforward 



Beginning with Equation 2.94, a sliding mode controller was formulated using 
standard SMC techniques, [Slotine 1991]. The sliding surface <j was defined as a 
function of the position error. 



<7 = 



r d_ 



■+■ X 



{ x ~x com )> 



(5.24) 



and the time derivative of o was defined as 

& = -T]sat(<j I <p) . (5.25) 

By defining the sliding surface in this manner, stability is guaranteed, based on Lyapunov 
analysis, since 

<j<7 < 0, Vf> 0. (5.26) 

Taking the time derivative of Equation 5.24 and equating it to Equation 5.25, the control 
input may be determined. 
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I i T 
n\n\ = - 





(5.27) 



T 



F + “ f -x com )-^ 2 (u r +u f -X c0m ) 



Using the signed square root of Equation 5.27, the commanded control input is found. A 
detailed description of this controller design approach may be found in [Riedel 1998b]. 

As seen in Equation 5.27, the commanded control input is a function of the 
system states, the fluid velocity (including the first and second derivative), the command 
inputs, and the "to-be computed" control input, due to the fact the system represented by 
Equation 2.94 is non-affine. To compute the required control input requires solving a 
difference equation in n, as well as measurements of the fluid velocity and its first and 
second derivative, making this control law extremely complex and possibly difficult to 
implement in real-time. To overcome these difficulties, some simplifications need to be 
made. 



If the thrust reduction term is ignored and treated as an unmodeled disturbance. 
Equation 5.27 reduces to. 



which requires only system states, fluid disturbance measurements and command inputs. 
To display how well this controller is capable of performing, again, consider the case of a 
monochromatic sine wave disturbance input, where the disturbance and its first and 
second derivative are known. When direct feedforward of the measured wave disturbance 
is embedded in the control system design, perfect cancellation of the wave disturbance 
effects on station keeping may be obtained. The simulated response of the PHOENIX, 
initially at five meters and closing to a commanded range of 0.5 meters, is displayed in 



a) Monochromatic Disturbance Input 




n\n\ = ■ 




(5.28) 
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Figure 5.24. Again, the results in Figure 5.24 are for demonstration purposes only, as a 
comparison with Case Ha. Perfect cancellation of the wave disturbances is not expected, 
since exact measurement of each wave disturbance component, i.e., u f ,u f and ii f , is not 



possible. The interesting point in this case is the propeller response. 
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Figure 5.24 Disturbance Cancellation Case Ilia, top to bottom respectively, position vs. 
time, propeller RPM vs. time, and a phase plane plot of the sliding surface 

When comparing the propeller response between the three cases that used a sine 
wave disturbance input, it can be seen that the SMC (Case Ilia) by far out performs the 
other designs. The position response is as desired, perfect cancellation, and the 
propulsion system is well within limits. This result is due to the fact that the system 
attempting to be controlled is highly nonlinear, requiring a nonlinear controller. 

b) PM Spectrum Based Disturbance Input 

Prior to continuing with any simulations to determine positioning 
performance, several simulations were conducted to determine the performance that 
could be obtained from the controller with and without all disturbance components 
available. Since the PM based disturbance input was generated using the techniques in 
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Chapter IV, the derivatives and phasing were known. Based on this, comparative 
simulations were conducted between controllers which used all the disturbance states and 
ones which used only the measurable fluid velocity state for disturbance rejection. 
Comparative results of one simulation are shown in Figure 5.25. 




Figure 5.25 Controller Performance Comparison, For A Controller That Uses All The 
Disturbance Components (Dashed Line), And A Controller That Uses Only Fluid 
Velocity For Disturbance Cancellation (Solid Line) 



As can be seen in Figure 5.25, the station-keeping improvements associated with 
including all components as opposed to including only the fluid velocity component is 
very small. In each case, the propulsion system response was within the vehicle's 
capability. As a result of the comparisons, it was determined that by using only the fluid 
velocity measurements, significant improvement with regard to positioning may be 
achieved. Based on this, the SMC takes the form 



n\n\ 



T 

~P 



6 - 2ccu r [au r \u r | + F]sign(u r ) + ... 

—+x com — 2X(au r \u r \ + F — ... 
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X C om)-^( U r+U f — X com ) 



(5.29) 
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which will be used for all remaining simulations. 

Using the PM based disturbance input allowed the SMC to be tuned to so 
that the controller would meet bandwidth requirements, limit propulsion system 
oscillations and avoid chattering. Controller parameters which provided a balanced 
design consisted of tj = 100, A = 1.0, and (f> = 0.5. The simulated position response of the 
Phoenix conducted with this SMC design is shown in Figure 5.26. 
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Figure 5.26 Position Response for Case Illb, PM Based Disturbance Input 

The position response shown in Figure 5.26 has a standard deviation of 
6.4 cm. This is twice as much as Case lib, Figure 5.18, however, the standard deviation 
of the propeller input for the SMC design is one-half a large as the LQR with disturbance 
estimator design, and is always within propulsion system limits. This can be seen in 
Figure 5.27. In addition, when comparing the two propeller responses, it appears that the 
SMC has a smoother output which will extend the life of the propulsion system. 
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c) Monterey Bay Disturbance Input 

Using the same design parameters that allowed the controller to achieve 
the performance depicted in Figure 5.27, the system was simulated with the disturbance 
input obtained from Monterey Bay. Since the disturbance magnitude of the Monterey 
Bay data is less than the PM based input, it is expected that the position response would 
also be less. By referring to Figure 5.28, it can be seen that this is in fact the case. 

The standard deviation of the position response has improved over the 
LQR based controller (Case lie) by a factor of 1.4 with only a slight increase in propeller 
rpms (5% compared to a maximum of 800). These results are shown in Figure 5.29. 

4. Disturbance Rejection Case Comparison 

After conducting the simulations for each of the cases with the various 
disturbance inputs, it was apparent that Case IE, the SMC with measured disturbance 
feedforward, out performed the other two cases and to most this is no surprise. What is 
interesting is the amount by which it out performed the other cases. 




Figure 5.27 Propeller Response for Case Illb, PM Based Disturbance Input 
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Figure 5.28 Position Response for Case IIIc, Monterey Bay Disturbance Input 
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Figure 5.29 Propeller Response for Case IIIc, Monterey Bay Disturbance Input 

Using the PM based disturbance input, simulations were conducted for each of the 
three case previously analyzed, LQR, LQR with disturbance estimator, and SMC with 
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disturbance measurement feedforward, with gains ranging from high to low. The attempt 
was to reproduce the "optimality" curve, Figure 5.15, for each controller to study the 
performance of each control solution. 

Conducting this study led to some very interesting results which are shown in 
Figure 5.30. As seen in this plot, the curves for each controller do not have the traditional 
optimal curve shape, i.e., as control input increases position error decreases. In fact, the 
curves indicate that for control designs ranging from low to medium gains, regardless of 
controller type, it would be better to have no control at all. 




Figure 5.30 Comparison Of Controllers For Various Gains, PM Based Disturbance Input 

The explanation for this can be seen in Figure 5.31. which superimposes 
the closed-loop vehicle frequency response, disturbance input to vehicle position output, 
over the disturbance spectrum for three different control gains, namely low medium and 
high. In Figures 5.30 and 5.31, point/curve "1" corresponds to a low gain, 
point/curve "2" to a medium gain solution and point/curve "3" is high gain control 

As Figure 5.31 displays, the low and medium gain solutions, with this 
particular disturbance, actually excites the vehicle, and not until a high gain solution is 
implemented does the vehicle actually reject the disturbance, Using this as an analysis 
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tool, the range of acceptable gains, for a particular disturbance input may be determined. 
In addition, it was quite evident that by feeding forward the measured disturbance using 
the SMC solution, that significant disturbance rejection was capable with much less 
power consumption. 
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Figure 5.3 1 Vehicle Frequency Response, (Disturbance Input To Position Output), 
Superimposed Over The PM Based Disturbance Input 



E. SUMMARY 

This chapter has outlined the various disturbance rejection techniques available to 
the control engineer. It has provided a summary of classical, modem and nonlinear 
control methodologies. Three case studies, which represent the basic design methods 
used to reject disturbances, were conducted and discussed for three different disturbance 
inputs. These studies showed that the SMC with measured disturbance feedforward is a 
far superior design approach for this particular class of problem, and provides significant 
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disturbance rejection performance for the same input power. Finally, an analysis 
approach that may be used to study gain selection and performance estimates has been 
presented. 
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VI. DISTURBANCE COMPENSATION CONTROLLER (DCC) 



A. INTRODUCTION 

This chapter will discuss the development of the real-time disturbance 
compensation controller (DCC) which will allow an AUV to dynamically position itself 
in the presence of waves. The chapter will begin with an overview of the DCC, followed 
by a discussion of an asynchronous Extended Kalman Filter for state and disturbance 
estimation. This nonlinear estimator is critical to the DCC performance since the SMC 
requires foil state feedback, and not all states are measurable. In addition, the EKF 
provides the controller with a smoothed estimate of the unmeasured fluid velocity which 
is used to compensate for the wave induced disturbance. 

Next, through the design and implementation of an asynchronous simulator, 
which realistically models the vehicle dynamics, the sensors including noise and the 
sensor processes, the DCC is tuned and the achievable performance is demonstrated. 

Lastly, it is shown that by properly weighting the noise covariance in the 
estimator the DCC reduces the transmission of sensor noise into the propulsion system 
while still maintaining the ability of the vehicle to hold position. 

B. DCC OVERVIEW 

The design of the disturbance compensation controller can be looked at as an 
optimization problem since there are competing goals. First, since the design 
requirement is to minimize position error in the presence of disturbances, a high gain 
control, as Chapter V discussed, is desirable. Using high gain control, the system 
becomes sensitive to measurement noise and uncertainty, thereby requiring the gain to be 
reduced to maintain stability. 

The estimator is needed to provide the unmeasurable states to the controller, and 
to filter the sensor noise thereby improving the systems performance. Here, the 
requirement is to accurately track the signal, again requiring a high filter gain, while 
smoothing the noise, (a low gain). As with the controller, trade-offs must be made. 
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The overall goal is to develop a combined controller/estimator which, when 
implemented, will enable the vehicle to maintain position while using noisy sensor 
information. The output of this system is a commanded voltage that is sent from the 
DCC process to the real-time execution computer, without excessive lags to ensure 
stability. A mathematical description to the above problem is given below, with a block 
diagram of the DCC in provided in Figure 6.1. 

State: x T - [X,u r , F\ d = u f 

System : x = f(x,n,d); y T =[x, ,]= Cx + DJ 

Disturbance: x f =Ax f +v u f =Cx f . (6.1) 

Control law : n = smc(x, u / ,x rom ) 

Estimator : \x,u f ]= EKF(f (x,n,d),A,y,n) 



Ground 

Velocity 




Figure 6. 1 Block Diagram of Disturbance Compensation Controller (DCC) 



C. STATE AND DISTURBANCE ESTIMATION 

There are many methods available to estimate states and disturbances in practice 
today. A few of these include the Luenberger Observer [Ogata 1990] and the Kalman 
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Filter [Gelb 1974] for linear systems, and the Sliding Mode Observer [Canudas De Wit 
1991], the "Rajamani" Observer [Rajamani 1998] and the Extended Kalman Filter for 
nonlinear systems. Each method has both pros and cons depending on the application. 
For this work, an Extended Kalman Filter was chosen since a relatively accurate vehicle 
model is available, and the stochastic nature of the disturbance. 

Kalman filtering is the process of recursively updating an estimate of systems 
states based upon measurements corrupted by noise. The system state is a collection of 
variables that describe the dynamics of a system, and in this case they are position, 
relative velocity and propeller thrust, of which only relative velocity is measurable. 

System states are updated with knowledge of system dynamics (vehicle model), 
measurement dynamics (measurement model), system noise (modeling uncertainty) and 
measurement noise (measurement errors). The system model is not perfect in describing 
the dynamics of the vehicle and will contain a certain amount of uncertainty, called 
system noise. There is also some uncertainty associated with each measurement taken. 
This uncertainty can be composed of both random white noise and a bias. Measurements 
which cannot be directly obtained, such as fluid velocity, are related to measurements 
which are directly obtainable, such as relative velocity and ground velocity, in the 
measurement model. Recursively updating means the Kalman filter does not need to 
keep record of all past measurements, only the most recent ones. 

1. Model and Filter Development 

Using the three state surge model developed in Chapter II, and a four state AR 
model for the wave dynamics, an augment state and disturbance model was formed, and 
used as the basis of an EKF. This model allows the disturbance to be treated as an 
additional state, where the vehicle states and disturbance estimates are filter outputs. The 
augmented vehicle and disturbance model is given by. 
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x = u r +u f 
u r = au r \u\ + F 




T T T 



T 




( 6 . 2 ) 




where the AR coefficients are found using the procedure outline in Chapter V. 

2. Kalman Filter Algorithm 

Using standard design techniques [Gelb 1974], the filter was developed and 
implemented using the following algorithm. First, the system model matrix A, system 
noise matrix Q, measurement matrix C, measurement noise matrix R, and the error 
covariance matrix P are initialized to appropriate values. The error covariance matrix can 
be thought of as a level of uncertainty in the state vector. Then the state vector, error 
covariance and measurement vector are propagated one time step using the model. 

When the new measurement is received, the innovation is calculated based on the 
difference between the measured values and the estimated values. Using the propagated 
error covariance, measurement noise matrix and measurement matrix, a gain is 
determined for the state vector and error covariance update. This process of propagating 
and updating is repeated through out the length of the vehicle mission. This recursive 
algorithm, in discrete form is given by, 







X k/k - 1 ~^k/k-i X k-\/k-l 




( 6 . 2 ) 
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where <I> represents the linearized system dynamics matrix, and h-C since the 
measurements are linear in the state. The continuous linearized matrices for this 
particular design are given as. 



A = 



C = 
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1 
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0 0 0 



(6.3) 



3. Asynchronous Data Processing 



In the preceding discussion, the data contained in the measurements was assumed 
to be received at the same time with equal intervals through out the mission. In reality, 
all measurements are not received at the same rate, therefore, the EKF design must allow 
for this asynchronous sampling rate. In the Phoenix AUV, the vehicle control loop 
currently runs at 8 Hz, while the RDI DVL runs at 2 Hz, and the SonTek ADV at 6 Hz. 
(See Appendix B for a more through description of sensor operations). The main data 
acquisition process samples the sensor processes at the same frequency as the control 
loop, however, if the sensor has not yet updated, the data acquisition process records the 
value of the previous time step. The filter allows for the varying measurement rates by 
using a dynamic switching of the measurement matrix, C, [Healey 1998]. The 
measurement matrix basically uses a zero-order hold on the measurement channel that 
has not been updated, and propagates the state using the previous measurement. 



D. ASYNCHRONOUS SIMULATOR DEVELOPMENT 

Using the filter design from the previous section, and the sliding mode controller 
developed in Chapter V, Equation 5.37, an asynchronous simulator was developed for 
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design validation. The simulator contains the non-linear vehicle dynamics, Equation 
2.94, asynchronous sensor models with measurement noise, seaway dynamics and the 
DCC. Using this simulator as a design tool allowed the DCC's control and estimation 
parameters to be adjusted prior to real-time implementation. Figure 6.2 shows a sample 
of the sensor outputs, during one of the simulation runs. As seen, the position output, 
which is a product of a navigation filter, is at 8 Hz. The relative velocity, u r , which is 
measured by the ADV, is at 6 Hz, and the ground velocity, u g , from the RDI is recorded 
at 2 Hz. In addition, the ADV output has noise imposed on the signal representative of 
the vendor advertised levels. 
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Figure 6.2 Asynchronous simulation with realistic noise models - Disturbance from PM 
Spectrum, H s = 1 meter, operating depth 10 meters 

Using this developed simulator, the DCC was adjusted to achieve an optimum 
design. The gains in both the controller and filter were adjusted so that performance 
requirements discussed earlier were met. Sample results showing the performance of 
final design are given in Figures 6.3 and 6.4. As can be seen, the estimates of both 
position and thrust track the actual values, and the position response is maintained within 
a standard deviation of 8 cm. This performance is extremely good, recalling that the 



130 



Longitudinal Position (m) 



same controller, with the same disturbance input was able to achieve a standard deviation 
in position of 6.5 cm, without noise and using full state feedback, see Figure 5.31. 




Time (s) 



Figure 6.3 Simulated and Estimated Position Response, Using Final DCC Design 




Figure 6.4 Simulated and Estimated Thrust Response, Using Final DCC Design 
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The stability performance of the estimator is shown through simulation, see 
Figure 6.5, since there are no formal proofs to determine the stability of combined 
nonlinear estimators and controllers. As seen, the error covariance levels all converge 
indicating a stable nonlinear filter design. Some of the covariance levels may appear to 
be "too high" giving the feeling that the filter is not properly designed, however, design 
decisions must be made to ensure that the filter lags are not too excessive, and that the 
estimator tracks well. 




Figure 6.5 DCC Error Covariance Evolution 



E. INITIAL IN-WATER TESTING 

Using short missions, that DCC was adjusted to achieve acceptable performance. 
These runs were performed on March 25, 1999, in Monterey Harbor. Of concern, was 
the amount of noise that was resident on the ADV sensor. This noise was far beyond the 
level which the vendor advertised. Using the design results from the simulations, the 
DCC was implemented in the Phoenix AUV. Figures 6.6-6.10 display initial results. As 
seen, the filter tracks the signals extremely well, including the noise. 
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Figure 6.6 Short Segment In- Water Results, Position for Radv=10 
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Figure 6.7 Short Segment In-Water Results, Relative Velocity for Radv=10 



133 



0.15 



0.1 



0.05 




3 



- 0.05 



- 0.1 



40 45 50 55 60 65 

Time (s) 





Figure 6.8 Short Segment In-Water Results, Fluid Velocity Estimate for Radv=10 
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Figure 6.9 Short Segment In-Water Results, Thrust Estimate for Radv=10 
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Time (s) 

Figure 6.10 Short Segment In-Water Results, Propeller RPMs for Radv=10 

This tracking of the noise has significant detrimental effects to the propulsion 
system as seen in Figure 6.10. The noise had been transmitted into the controller 
resulting in severe oscillation in the propeller response. These oscillations eventually 
lead to mechanical failure of the propulsion system shafting due to the shearing of 
connecting pins. 

Using the information obtain during this set of runs allowed the filter gains to be 
adjusted to eliminate the transmission of sensor noise into the controller. Using linear 
design techniques, the combined controller filter transfer function from ADV input to 
propeller output was formed. By adjusting the level of the measurement noise 
parameters, attenuation of the noise into the control system was accomplished. These 
results are shown in Figures 6.11 and 6.12. As the Figures shown, for the smaller noise 
estimate (Radv=10), the noise transmission is amplified over most of the range of the 
controller, while for the higher noise estimate (Radv=100), the input to the controller in 
attenuated. This improvement in frequency response will reduce the propeller 
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oscillations, thereby minimizing the chance of mechanical failure of the propulsion 
system. 
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Figure 6.1 1 DCC Frequency Response, ADV Input to Propeller Output , Radv=10 
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Figure 6.12 DCC Frequency Response, ADV Input to Propeller Output , Radv=100 
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Using the new design values, the DCC was again tested in Monterey Harbor. The 
results of this testing are shown in Figures 6.13-6.17. 
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Figure 6.13 Short Segment In-Water Results, Position for Radv=100 
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Figure 6. 14 Short Segment In-Water Results, Relative Velocity for Radv=100 
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Figure 6.15 Short Segment In-Water Results, Fluid Velocity Estimate for Radv=100 
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Figure 6. 16 Short Segment In- Water Results, Thrust Estimate for Radv=100 
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Figure 6.17 Short Segment In-Water Results, Propeller RPMs for Radv=100 

As the result of the tuning of the DCC, the performance has improved 
dramatically. As before, the DCC maintains position extremely well, with a much 
reduced propeller response. Comparing the magnitude of the estimated fluid velocities 
between the two designs. Figures 6.8 and 6.15, it can be seen that for the same magnitude 
of input disturbance, position response has remained unchanged, but propeller response 
has reduced increasing the life of the propulsion system. 

F. SUMMARY 

The design and implementation of a new Disturbance Compensation Controller 
(DCC) has been presented. The results indicate that by using a properly tuned system, 
the ability of a tetherless underwater vehicle to dynamically position itself with minimal 
input is possible. Although no formal proof for stability is available, asynchronous 
simulations have demonstrated that the DCC is stable and provides acceptable tracking 
and estimation of state and disturbance inputs. 
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VII. ESTIMATION OF WAVE DIRECTIONALITY FROM A MOVING 

PLATFORM 



A. INTRODUCTION 

This chapter will outline the underlying principles used in identification of wave 
directions from standard wave following buoys. It will present the mathematical 
formulas used in determining the wave direction as a function of frequency. Extension of 
these algorithms to subsurface velocity sensors will be made, where, through the use of 
the Doppler equation, a moving AUV can be used to determine wave directions. Lastly, 
it will be shown how a control command can be obtained from the frequency dependent 
wave direction estimates. 

The information in this chapter is not new, only the application to which this 
method is applied. For more detailed descriptions of the mathematical formulations 
presented in this chapter, the reader is referred to [O'Reilly 1996} and the references 
therein. 



B. WAVE SPECTRA AND DIRECTIONAL ESTIMATES 

As discussed in Chapter II, a wave record, Tj (t), measured at a fixed location can 
be represented as a linear superposition of waves of different frequencies. Wave 
components with different frequencies are usually assumed to be statistically independent 
because they are generated by random wind forces at different locations. From the 
central limit theorem it follows that the probability distribution of T](t) is approximately 
Gaussian, consistent with many observations, [Soong 1993]. 

The procedure presently employed by many of the operational data buoys in 
based on Fourier analysis. In Fourier analysis it is convenient to work with complex 
exponentials rather than sine and cosine functions, there fore using the relation 

cos(«-h») = «P«<a+0))+«xp (-<«*+»)) (7.J) 

the expression for the surface elevation can be written as 

7j(t) = ^A 0) exp (icot) , (7.2) 
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where 



K = ^exp(^J 

and the summation is over both positive and negative frequencies. 

As discussed in Chapter IE, the energy spectrum E( (O), is defined as 



E(co) = 



A(0 



(7.3) 



(7.4) 



where ( ) indicates an average over many data records and Aco is the spacing of the 



frequency bands. The spectrum is closely related to the energy of the waves, and 
represents the distribution of wave energy as a function of frequency. 

To describe the spatial and temporal characteristics of the sea surface linear 
superposition of wave components is used. In exponential terms this can be represented 
as 



J](x, y, t) = 2 X A o>.e exp(i(£(*cos 0 + y sin 6) - wt )) (7.5) 



0 ) 6 



where x, y are the horizontal spatial coordinates, and (0 and k obey the dispersion relation. 
The frequency directional wave spectrum is defined as 

|2 ' 



E«o) = 



co 7 0 



A(0 Ad 

and describes the distribution of energy as a function of frequency and direction. 



(7.6) 



C. WAVE BUOYS 

The most commonly used instrument for measuring waves in deep water is the 
"heave, pitch and roll buoy" that measures the surface height and slope in two orthogonal 
directions. The newer Datawell® Directional Waverider measures 3-component 
accelerations of the buoy which are integrated to yield the horizontal and vertical 
displacements of the buoy. The hull and mooring were designed to give the buoy good 
wave following characteristics, thereby allowing the buoy displacements to approximate 
the displacements of an actual water particle at the sea surface. 
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For a wave train propagating in the positive ^-direction, the fluid particle motion 
is given by Equation 2.49. For the more general case of a wave train propagating at some 
angle relative to the x-axis, it can be shown that the flow field is given by 

u(x, y,t) = ao)cos0txp(kZ)cos(k(xcosd+ y sin d)-ojt) 

v(x, y,t ) = acosvn #exp(£Z)cos(£(xcos# + y sin 6) -OX) (7.7) 

w(x, y,t) = a<yexp(£Z)sin(fc(xcos# + y sin 6) — OX) 

Let the average position of the buoy be given by x=y=z=0. For small amplitude waves, 
the expected buoy displacements are small compared to the surface wavelength, therefore 
the buoy motion can be approximated by the fluid velocity at x=y=z= 0. 

For a full spectrum of waves, the buoy displacements can be expressed using 
complex notation as 

— iA ae cos#exp(in?) 

o) e 

m= IS - iA m6 sin 6 exp(/M) (7.8) 

co 6 

Z (*) = Z Z “ A co,8 exp(iwt ) 

c o B 

where the -i is due to the 90° phase difference between the vertical and horizontal 
displacements. The expressions in Equation 7.8 can be written using Fourier transforms 
as 

X (0 = ^ X ( a )) exp(iwt) 

0) 

Y(t) = Y J Y(0))txp(iwt) , (7.9) 

CO 

Z(t) = ^Z(ty) exp(/wr) 

CO 

where the Fourier transforms are given by 

X (°>) = 'Z<- iA co,e cos # 

6 

Y (<*>) = ^-iA^e sin6 . (7.10) 

6 

Z{(D) = YjKe 

8 

To derive the relationships between the measured time series and the unknown 
frequency-directional wave spectrum the cross spectrum must be considered. In general, 
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the cross spectrum between two time series X(t) and Y(t) with Fourier transforms X(co) 
and Y( co) is defined as 



c yy (ty) — 



( X(6))Y'(0 ) )) 



Act) 



(7.11) 



where * indicates the complex conjugate, [Soong 1993]. Substitution of Equation 7.10 
into Equation 7.11 yields 



C XY (co) = ^ cos ^ sin 6E{co,6) (7.12) 

e 

where it is assumed that the wave components propagating in different directions are 
statistically independent. The cross spectrum Cxy can be expressed in continuous form as 



27! 

C XY {CO) = J cos 6 sin 6E{co , 6)d6 . (7. 1 3) 

o 

Cross-spectra of the other time series pairs can be obtained in a similar manner. The full 
set of relations for buoy cross-spectra can be found in [Dean 1984]. 

It is convenient to define a normalized directional distribution of energy at a 
frequency co as 



S{6\co) 



Ejo^d) 
E{t o) 



(7.14) 



with unit integral 



In 

\S{6-co)dd = 



In 

\E{co,d)dd 



E{co) 



= 1 . 



(7.15) 



E{co) E{co) 

With this definition. Equation 7.13 and the other referenced spectral relations can be 
combined and expressed in terms of S{6\co) . Dropping the frequency dependence these 
relations can be expresses as 
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Im (C xz ) 



((Cn+Cyy'Pn )" 2 

Im (Cyz ) 



((Cxx Cyy )C a ) 



1/2 



= Jcos 6S(0)dd = a y 
o 

2n 

= Jsin 6S(d)dd = b l 



C -C 
^ xx ^ 



0 

In 



YY 



Cxx + CyY 
2 Re(C m ) 
Cxx + Cyy 



= jcos26S(0)d$ =a 2 

0 

In 

= | sin 29S(9)dO = b 2 



(7.16) 



These four relations between the cross-spectra of the buoy measurements and the 
directional distribution of wave energy, derived by [Long 1980], form the basis for most 
of the buoy analysis techniques currently used. 



D. EXTENSION TO SUBSURFACE SENSORS 

As discussed in the previous section, the motion of a wave buoy is directly related 
to the fluid velocity, therefore, the cross-spectra of a tri-directional current meter yields 
the same low resolution directional wave information obtained from buoy measurements. 
Substituting the normalized spectra of the vertical (Z) velocity component w, and the 
horizontal ( X , Y) velocity components u and v into Equation 7.16, the lowest four Fourier 
moments of the directional distribution of wave energy can be obtained and are given by 



1 [Create uu (a)) + C n (a))Y ’ 


(7.17) 


b(m')- Im (C m (co)) 

1 [c^icotc^ + c^T' 


(7.18) 


C (o))-C (CO) 
aJcd)= uuK ’ nK } , 
C„(a>) + C„(a) 


(7.19) 


2Re(C » <ft,)) , 
C uu (co)+C n (cd) 


(7.20) 



where C( on) is the spectral matrix of the velocity components u, v, w. Since the wave 
direction, 6 , is referenced to the navigation frame (N-E-D), vehicle borne sensor 
measurements must be transformed prior to use. It is interesting to note that the estimates 
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of these directional moments are insensitive to errors, so long as the errors are the same 
on all measurement axes of the sensors, which is typical with oceanographic sensors 
installed on AUVs. 

The objective of the data analysis is to infer the directional distribution S( 6), from 
the four measured moments a 2 , bi, a 2 , and b 2 . The most widely used techniques are 
described below. 

a) The Cosine Power Distribution 

[Longuet-Higgins 1963] introduced a simple cosine-power distribution, 



with Omean the mean propagation direction, 5 a parameter that controls the width of the 
distribution and A, a normalization coefficient. The parameters 6 mean and s are 
determined by fitting Equation 7.21 to the relations given in Equations 7.17-7.20. 

The main drawback to this simple method is that Equation 7.21, with only two 
free parameters can describe only unimodal distributions, and thus fails in situations with 
a bimodal sea state (e.g., multi-directional seas during a wind veering event or swell 
arriving from two different sources). 

b) The Maximum Entropy Method 

[Lygre and Krogstad 1986] introduced the maximum entropy or MEM 
estimate of S(9). Unlike Equation 7.21, this approach can describe both unimodal and 
bimodal distributions and exactly fits the relations given in Equations 7.17-7.20. This 
directional spectrum is given by 



S{9) = Acos 2s 




(7.21) 



S(9) = 



1 1 01^1 02^2 



(7.22) 




with 
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(7.23) 



Cj = a x + ib x 
c 2 = a 2 +ib 2 



<p =£iZ££l 

i-hl 

02 = C 2 ~ C \ < P\ 



Still, the directional distribution is poorly constrained by only four moments and the 
estimates require careful interpretation, [Krogstad 1991]. 



c) Mean Direction and Directional Spread 

An alternative approach that avoids the pitfalls of S( 6) estimation, is to 
describe the directionality of waves by a few simple parameters. For narrow S(0), a 
mean propagation direction 6 m and a root-mean-square measure of the directional 
spreading energy a e can be defined in terms of the first-order and second-order Fourier 
moments a,, b x , a 2 and b 2 [Kuik 1988], given by. 



6_ = tan 



-i 



V 

y a \ j 



a] = 2[l - [a, cos(0 m ) + b x sin(0 m )J 



6 m = — tan 1 

m 2 



(u \ 



\ a 2 J 



(7.21) 

(7.22) 

(7.23) 



<j] = \b -\- a 2 cos (26 m )+b 2 sin(2^ m )U. (7.24) 

Again, this method fails to identify bimodal distributions but it is useful to 
determine a base direction so that a control command could be determined. More on this 
approach will be discussed later in this chapter. 



E. CORRECTION FOR A MOVING PLATFORM (THE DOPPLER 
EQUATION) 

The equations for the wave directionality estimations presented in the previous 
section is valid for a non-moving sensor. However, when information is obtained from 
an AUV, corrections must be made to account for the vehicle motion. As discussed in 
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Chapter III, the wave frequency which the vehicle encounters while moving through the 
wave field has been shifted. This shift can be determined by using the well known 
Doppler equation. Equations 2.56 or 4.25. 

Using the techniques outlined in the previous section will give the wave 
directional distribution as a function of vehicle encounter frequency. If these estimations 
are used in conjunction with the Doppler equation in a recursive manner, the estimation 
of S( 6) can be found as a function of true frequency. 

The method by which this is performed is outlined below; 

• Determine the three components of the fluid velocity in vehicle body fixed 
coordinates. 

• Transform the fluid velocities into the global navigation frame using Equation 

2 . 8 . 

• Compute the auto- and cross-spectra of the fluid velocity components. 

• Determine the Fourier moments using Equations 7.17-7.20. 

• Convert the Fourier moments into Krogstad notation and use the MEM to 
determine the directional distribution S( 6). 

• Using the vehicles mean velocity, and the mean direction obtained from 
Equation 7.21 or 7.23, apply the Doppler equation to determine the frequency 
shift. 

• Return to 3 and complete the process until frequencies converge. 

The corrections to the estimation of S(6) using the Doppler equation are quite 
small for slow moving vehicles. Considering, for example, the NPS Phoenix AUV which 
has a maximum forward velocity of 1.5 m/s, the error associated with not using the 
Doppler equation are approximately ± 1 sec in identification of wave periods between 4 
and 40 seconds. Similarly, the error in direction estimation is approximately 5-7 degrees. 
When an AUV goes into a station keeping mode, where vehicle velocities are 
significantly reduced, the modifications required due to the Doppler shift are negligible. 
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F. 



DETERMINATION OF CONTROL COMMANDS 



The resulting directions and variances, from Equations 7.21-7.24, are for each 
frequency component of the wave field. To use this information in determining a heading 
command, single direction must be found. Bulk Fourier moments, weighted by the 



energy density of the wave field, 



'*r 




'a^co) 




1 fu 
C /; 


b,{(0) 

a 2 (co) 


M. 




_b 2 {co)_ 



dco 



with E b the swell variance given by. 



/. 

E‘ ' = f E(a>)da> , 



(7.2#) 



(7.2#) 



may be substituted into Equations 7.21-7.24 to yield a bulk fluid direction and variance. 
It is this bulk fluid direction that is used in generation the heading command dependent 
on the mission requirements. 



G. SUMMARY 

This chapter has presented the techniques currently employed for the 
determination of wave directional estimations from standard wave following buoys. It 
has extended this analysis for use in determining directional estimates from data gathered 
by an Autonomous Underwater Vehicle. Using this data gathered, an approach was 
presented which will allow the deployed vehicle to obtain information about the 
directionality of its working environment thereby allowing the vehicle to have 
information available to make decisions regarding mission execution. 
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VIII. EXPERIMENTAL RESULTS 



A. INTRODUCTION 

With the theoretical presentation of this dissertation complete, the primary focus 
of this chapter will be on the experimental validation of the Disturbance Compensation 
Controller (DCC). In addition, results from the ONR sponsored AUVFEST ’98 will 
prove the concept of wave direction estimation from a moving platform presented in 
Chapter VII. 

This chapter will begin with a discussion of the real-time implementation of the 
DCC in the NPS Phoenix AUV. It will follow with a presentation of the experimental 
results, conducted in Monterey Harbor, which displays the achievable performance of the 
DCC. Comparison between theoretical and experimental results will be made. Lastly, 
directional energy and spectral plots obtained by Phoenix, while operating in the Gulf of 
Mexico, south of Gulfport, MS, will be shown. 

B. SOFTWARE IMPLEMENTATION 

The implementation of this control process is unique since it is split between the 
two CPUs installed in Phoenix. The NPS AUV uses a Pentium based PC- 104 running 
QNX and a GESPAC Card Cage running OS9 for mission control and execution The 
DCC requires information from both processors, connected by Ethernet sockets, to 
compute and pass the commanded propeller RPMs to the execution level. 

The control architecture presently running in Phoenix is based on shared memory 
processes. The PC- 104 computer runs a “main” process that controls the synchronization 
of the data sharing, while the GESPAC clock controls the real-time control features. The 
two-processors use the shared memory as the common data buffer, controlled by 
semaphores to ensure the information transfer is consistent with the clock speed. A 
graphical representation of this description is shown in Figure 8.1. As seen in the 
graphic, for the DCC implementation, all needed processes are run in the PC- 104 with the 
main purpose of the GESPAC is to provide the commanded voltages to the propulsion 
motors. 
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Figure 8. 1 Software Implementation of DCC 

A block diagram of the DCC implementation in the Phoenix AUV was given in 
Chapter VI, Figure 6.1. It is redisplayed as Figure 8.1 for easy reference. This diagram 
represents the melding of the software and the hardware in the vehicle. The ground 
velocity is from the RDI, the relative velocity from the ADV and y/, r from the 
directional gyro. 



Ground 

Velocity 




Figure 8.2 Block Diagram of DCC Implementation 
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c. 



EXPERIMENTAL VALIDATION OF THE DCC 



The DCC was tested in Monterey Harbor between the months of March and May 
1999. During this time, the Phoenix was held under surge control for over 90 minutes, 
during various runs, without a drive off. Table 8.1 provides a sample of the runs 
conducted during the validation of the controller. 



Date 


Run # 


Length 


DRK 


Oil /it max 


comments 


4/2/99 


4 


4 min 


0.3624 


2.96 


high gain, 
short run 




5 


4 min 


0.6324 


3.08 


high gain, 
short run, 
vehicle physical 
disturbed 




6 


4 min 


0.4312 


0.265 


high gain, 
short run 




8 


4 min 


0.5090 


0.285 


high gain, 
short run 


4/22/99 


3 


10 min 


0.5508 


0.108 


high gain, single 
shaft 


5/25/99 


6 


10 min 


0.3620 


0.192 


medium-high gain, 
ADV noise 
problem 




8 


10 min 


0.3978 


0.126 


medium-low gain, 
ADV noise 
problem 




9 


10 min 


0.4957 


0.083 


low gain, ADV 
noise problem 




11 


10 min 


0.3587 


0.202 


medium-high gain, 
ADV noise 
problem 




12 


10 min 


0.4276 


0.144 


medium-low gain, 
ADV noise 
problem 



Table 8.1 Sample Summary of DCC Validation Runs 
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Defining a measure of performance, the disturbance rejection ratio ( DRR ), as the 
ratio of standard deviation of the vehicle ground velocity to the standard deviation of the 
fluid velocity, the ability of an AUV to reject disturbances for different conditions and 
control designs can be compared. Referring to the DRR definition, for perfect 
disturbance cancellation the DRR will be equal to zero, while for designs where the 
control input excites the vehicle, as was shown in Chapter V, the DRR will be greater 
than one. For each operating point, the standard deviation of the propeller response is 
normalized by the maximum propeller revolutions, 

Table 8. 1 indicates that excellent disturbance rejection was achieved, even for the 
short runs where only limited statistical information was recorded. The tests showed that 
even when the vehicle was disturbed by a source other than the fluid velocity, it was able 
to return to the commanded position in a stable fashion. 

A series a plots. Figures 8. 3-8. 8, show the results of one of the validation runs. 
This run was conducted on April 22, 1999 in Monterey Harbor. The Phoenix was 
commanded to a navigational position of 0 meters in the longitudinal direction. As the 
results indicate, the vehicle behaved as expected. The standard deviation of the 
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Figure 8.3 Comparison of Measured and Estimated Position, April 22, 1999, Run # 3 
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Figure 8.4 Propeller Response, April 22, 1999, Run # 3 
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Figure 8.5 Measured Ground Velocity, April 22, 1999, Run # 3 
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Figure 8.7 Fluid Velocity Estimate, April 22, 1999, Run # 3 
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Figure 8.8 Estimated Thrust, April 22, 1999, Run # 3 



positional error was 9.6 cm with ground velocity standard deviation of 1.5 cm/s. 

This run was the most interesting of the validation runs conducted. At the 
beginning of this run, it was noticed that the starboard shaft was not turning. Even with 
this propulsion system casualty, the vehicle was able to hold position and the controller 
did not go unstable. With only one shaft turning the effective input gain for the vehicle 
was reduced by 50%. Operations of this nature indicate a very robust design. It can be 
seen in Figure 8.4, that there is a small increase in propeller revolutions around the 50 
second point of the run. Data analysis indicated that this was approximately when the 
starboard shaft failed. Investigation into the cause of the shaft failure determined that a 
universal joint in that shafting had worked loose. 

As a graphical representation of the performance expected for the DCC a various 
simulations ware conducted, using the asynchronous simulator discussed in the Chapter 
VI, with the estimated fluid velocity obtained during this run (April 22, 1999, run # 3) as 
the disturbance. The gains on the DCC were varied to produce a position response verses 
propeller response curve. The actual experimental results, presented in Table 8.1, were 
superimposed on the theoretical curve obtained from simulation These results are shown 
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in Figure 8.9. As seen, the experimental and theoretical results are very close indicating a 
physically realistic simulator. 

The comparisons displayed in Figure 8.9, contrast data analyzed from 
experimental results together with computer predicted results with the same disturbance, 
and yield insight into the achievable performance of the DCC. It indicates that there is a 
limit to the amount of disturbance rejection that is physically realizable. This limit is 
controlled by ability of the propulsion system to produce the needed input to maintain 
position without excessive oscillations. The excessive oscillations have a detrimental 
effect of the life of the propulsion system. 

As a note, the short runs, displayed in Figure 8.9, were conducted with a 
controller gain parameter of rj = 100, a high gain. If the length of these runs were 
extended, these points would move closer to the curve as with the other runs displayed. 



-to 



c 

o 

•a 




Figure 8.9 Comparison of DCC Performance, Simulation and Experimental 

Up to this point, the only results presented are for the Phoenix maintaining 
position to the origin, the point which the run was initiated. Questions arise as to how 
effective the controller is in dealing with transients. This question may be answered by 
referring to Figure 8.10. This figure depicts the transient response of the Phoenix for the 
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various DCC gains presented in Figure 8.9. As the figure indicates, the controller deals 
with transients extremely well. 

The responses displayed in Figure 8.10 are for the regulator solution. What is 
meant by this, recalling that the SMC formulation requires kinematically consistent 
position, velocity and acceleration, is that no command inputs, other that position were 
used. In doing this, it is expected that the vehicle will over shoot and oscillate around the 
commanded position consistent with some settling time. 
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Figure 8. 10 Comparison of Transient Response for Various Control Gains 

With these transient responses come some issues with regard to operational 
implementation. If the goal is to position the vehicle close to, but with out touching, an 
object, some means of predicting the transient must be available. A resulting question 
that needs to be answered is; Does the developed simulator, which, based on the 
comparison in Figure 8.9, accurately predict the transient response? By comparing the 
results of the experimental runs and the simulated results, for the same disturbance input 
and DCC design, see Figure 8 11, the question can be answered, "yes.” 

As seen in this plot, the simulated results accurately reflects the measured 
transient response of the Phoenix. The steady state response, however, does not match 
exactly. The reason for this is two-fold First, the Phoenix, for recovery reasons, is 
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maintained approximately two-pounds buoyant. This weight and buoyancy mismatch, as 
discussed in Chapter II, cause additional excitation forces resulting from the wave 
induced fluid accelerations. Since the fluid acceleration cannot be measured, this 
additional excitation force is difficult to replicate in simulation yielding errors between 
the real and simulated response. Second, the experimental results are measured from a 
6DOF rigid body, where as the simulator results come from a 1DOF surge model. The 
coupling effects from the surge-pitch dynamics will effect the comparison. 




Figure 8. 1 1 Transient Response Prediction of the DCC 



D. WAVE DIRECTIONAL ESTIMATION USING THE NPS PHOENIX AUV 

During November 1998, the NPS Center for AUV Research, under the direction 
of Professor Anthony Healey, participated in the ONR sponsored AUVFEST '98. This 
AUV technology demonstrator was held in the Gulf of Mexico, south of Gulfport, MS. A 
complete description of the exercise can be found in Appendix D. 

The Phoenix AUV was used during this exercise as a mobile sensor to gather 
oceanographic data. Using the concepts presented in Chapter VII, the vehicle conducted 
27 short-term sampling missions. The products that were obtained included directional 
energy plots, directional spectrum plots and mean current estimations. 



160 



The key to the ability of Phoenix to obtain data capable of producing this 
information is the combined ADV/RDI/MotionPak/Directional Gyro sensor suite. With 
these sensors, accurate three dimensional fluid velocities, expressed in global quantities, 
were capable of being obtained in post-processing. Since the RDI/ADV sensors are 
collocated, little vehicle induced motion remains after processing the data. 

The data obtained validated the concept of obtaining tactical oceanographic data 
from an underwater vehicle. During the collection of the data, remnants of Hurricane 
Mitch were still present in the Gulf, providing an excellent source of ground swell. In 
addition, there was a significant wind generated wave component in a different direction 
than the swell component, resulting in a multi-modal spectrum. 

The results presented in Figures 8.12-8.14 provide a sample of the oceanographic data 
obtained during this offshore exercise. As seen in Figure 8.12, the bi-modal properties of 
the seaway are captured, as well as an estimate of the significant wave height. The ability 
to estimate the dual directions is due to the use of the MEM algorithm presented in 
Chapter VII. Figure 8.13 presents the associated spectrum plots for this energy estimate. 
The ability of an AUV to estimate currents is shown in Figure 8.14. Using a triangular, 
time based run, the current can be determined using set and drift calculations from the 
error in final position as well as the heading error on each leg. This information can be 
feed directly into the vehicle’s navigation process to account for the offset due to current 
thereby increasing navigation accuracy. Short term averages from each of the three legs 
are in general agreement with the overall average determined from the navigational drift. 

E. SUMMARY 

This chapter has presented the experimental results associated with this 
dissertation. The results validate the theory and show that it is possible for a small AUV 
to reject wave induced disturbances making it capable of performing positioning tasks in 
shallow water. In addition, the robustness of the design to sensor noise and propulsion 
faults has been demonstrated. 
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Phoenix Survey Data 



Significant Wave Height (m): 0.351 
Peak Frequency(Hz)/Period(sec): 0 25 / 4 
Peak Direction: 90degrees 

Max. Energy: 0.005203cm 2 /Hz 
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Figure 8.12 Sample Direction Energy Plot From Phoenix Data, Nov. 4, 1998 (Run # 9), 

Gulf of Mexico 
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Figure 8. 13 Sample Direction Spectrum Plots From Phoenix Data, Nov. 4, 1998 

(Run # 9), Gulf of Mexico 
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As a supplement, due to the sensor suite available, it was shown that tactical 
oceanographic data is obtainable from a moving AUV. The vehicle in this manner 
becomes an intelligent, mobile off-board sensor, thus presenting the argument for AUV 
deployment with operational fleet units. 




-200 -150 -100 -50 0 50 100 150 200 
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Figure 8.14 Short-term Current Estimation from Phoenix, November 8, 1999, (Run # 2), 

Gulf of Mexico 
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IX. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

This research has shown through modeling, simulation and experimental 
validation that intervention tasks performed by intelligent underwater robots are 
improved by their ability to gather, learn and use information about their working 
environment. The development, implementation and verification of a real-time 
embedded Disturbance Compensation Controller (DCC) for small AUVs, has provided a 
new technology showing that it is possible to use underwater vehicles for station-keeping 
tasks in shallow water. 

The work conducted in support of the dissertation objectives has provided three 
specific contributions to the field of tetherless underwater robotics. First, a new 
generalized approach to the modeling of small underwater vehicles subject to shallow 
water wave and current effects was developed. Using appropriate modeling 

formulations, the fluid disturbance effects are incorporated directly into the equations of 
motion leading to model based control laws for disturbance compensation. In addition, 
this formulation provides the ability to construct a generalized distributed simulation 
capability for the evaluation of underwater vehicle systems in shallow water, which is 
particularly useful to U.S. Navy tactical decision making. 

Secondly, it is proven using asynchronous simulations and in-ocean experimental 
validations, that substantial compensation of wave induced disturbances may be achieved 
from direct on-line measurements of the water column velocities. This technique 
eliminates the need to develop and incorporate sophisticated predictive disturbance 
models in the control system design. 

Thirdly, it is shown how small underwater vehicles, using direct fluid 
measurements, can obtain short-term wave magnitude and direction, as well as current 
estimates, thereby providing useful information in the area of tactical oceanography. It is 
also shown how a general seaway direction may be obtained from this information for 
use in mission planning and control. 
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In addition to the general contributions listed above, several specific conclusions 
may be drawn from this research. In particular, 

• The input requirements associated with the DCC are vehicle position, relative 
velocity, propulsion force, and water column velocity. The attractive nature 
of the DCC is that these quantities may be measured or estimated from a hull 
mounted sensor suite. With these the full benefits of DCC can be realized. 

• The DCC proved to be very robust to sensor noise and propulsion system 
faults. Stable vehicle performance was maintained in the presence of the loss 
of a propulsion shaft. 

• The design of the propulsion system must allow for a rapid, oscillatory 
response to the call for propeller speed commands. In addition, considerations 
for fault tolerant operations are recommended. 

• For the application of dynamic station-keeping, propulsion system lags and 
associated thrust reduction terms must be identified and taken into account. 

• Identified parameters of the nonlinear model provided stable and easily 
tunable controllers, as verified by in-water experiments. Excellent agreement 
was achieved between experimental and simulated results. 

• Seaway models developed using AR representations require high order to 
effectively achieve spectral matching. However, lower order models (four or 
six) can be used adequately for estimation and control. 

• Extended Kalman filtering methods for seaway estimation appear to be 
satisfactory. 

B. RECOMMENDATIONS FOR FUTURE RESEARCH 

As a result of the work performed in this dissertation several research areas have 
arisen requiring further investigation. These include: 

• The validation of the 6DOF model for other control modes with field data is 
recommended. This is particularly interesting in flight control where motion 
minimization is critical to improve side-scan imagery. 

• Application of the disturbance compensation techniques presented herein to 
other control modes used in shallow water positioning needs inquiry. By 
achieving compensation in all three planes, small AUVs may be employed in 
a number of positioning tasks, including mine neutralization. 
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• A sensitivity study as to the required degree and accuracy of the disturbance 
model, with regard to disturbance compensation, is warranted. It is believed 
that improved disturbance dynamics models may increase DCC robustness. 

• While this work has used EKF methods for the identification of seaway 
dynamics, other techniques, such as neural networks and MUSIC, may have 
advantages that could be explored. 

• The DCC formulation does allow for prediction of future water column 
velocities. This added information may possibly provide additional benefits 
not explored here, although, this study has indicated that the DCC is not 
highly sensitive to future information. 
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APPENDIX A. EQUATIONS OF MOTION AND PARAMETERS FOR THE NPS 

PHOENIX 

The equations of motion, parameter description and parameter values used to 
simulate the dynamic behavior of the NPS Phoenix AUV is given in this Appendix. 

PHYSICAL PARAMETERS 



Parameter 


Description 


Value 


W 


Weight 


1934.9 N 
(435 lbs) 


m 


mass 


197.2 kg 

(13.52 lb-sec 2 - ft' 1 ) 


B 


Buoyancy 


1934.9 N 
(435 lbs) 


l 


Characteristic Length 


2.225 m 
(7.302 ft) 


Ixx 


Mass Moment of Inertia 
about x-axis 


3.66 N-m-sec 2 
(2.7 ft-lb-sec 2 ) 


lyy 


Mass Moment of Inertia 
about y-axis 


56.94 N-m-sec 2 
(42. ft-lb-sec 2 ) 


I zz 


Mass Moment of Inertia 
about z-axis 


61.01 N-m-sec 2 
(42.0 ft-lb-sec 2 ) 




Cross Product of Inertia 
about xy-axes 


0.0 N-m-sec 2 
(0.0 ft-lb-sec 2 ) 


hz 


Cross Product of Inertia 
about yz-axes 


0.0 N-m-sec 2 
(0.0 ft-lb-sec 2 ) 


hz 


Cross Product of Inertia 
about xz-axes 


0.0 N-m-sec 2 
(0.0 ft-lb-sec 2 ) 
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XG 


x Coordinate of CG From 


0.003 m 




Body Fixed Origin 


(0.01 ft) 


yc 


y Coordinate of CG From 


0.003 m 




Body Fixed Origin 


(0.01 ft) 


Zc 


Z Coordinate of CG From 


0.0128 m 




Body Fixed Origin 


(0.042 ft) 


xb 


x Coordinate of CB From 


0.003 m 




Body Fixed Origin 


(0.01 ft) 


ys 


y Coordinate of CB From 


0.00 m 




Body Fixed Origin 


(0.0 ft) 


Zb 


z Coordinate of CB From 


0.00 m 




Body Fixed Origin 


(0.0 ft) 


Xbvt 


Location of Bow Vertical 


0.432 m 




Thruster from CG 


(1.420 ft) 


%svt 


Location of Stem Vertical 


-0.432 m 




Thruster from CG 


(-1.420 ft) 


Xbit 


Location of Bow Lateral 


0.585 m 




Thruster from CG 


(1.920 ft) 


Xslt 


Location of Stem Lateral 


-0.585 m 




Thruster from CG 


(-1.920 ft) 


yis 


Location of Left Propeller 


-0.10 m 




from CG 


(-0.33 ft) 


yrs 


Location of Left Propeller 


0.10 m 




from CG 


(0.33 ft) 
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CONTROL INPUTS 



Parameter 


Description 


Value 


Fis 


Left Propeller Force 


±44.45 N 






(± 10 lbs) 


Frs 


Right Propeller Force 


±44.45 N 






(± 10 lbs) 


Fbu 


Bow Lateral Thruster 


± 22.25 N 




Force 


(±5 lbs) 


F s it 


Stem Lateral Thruster 


± 22.25 N 




Force 


(± 5 lbs) 


Fbvt 


Bow Vertical Thruster 


± 22.25 N 




Force 


(± 5 lbs) 


Fsvt 


Bow Vertical Thruster 


± 22.25 N 




Force 


(± 5 lbs) 


Sbr 


Bow Rudder Deflection 


± 0.4 rad 


dsr 


Stem Rudder Deflection 


± 0.4 rad 


Sb P 


Bow Plane Deflection 


± 0.4 rad 


Ssp 


Stem Plane Deflection 


± 0.4 rad 
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NON-DIMENSIONAL HYDRODYNAMIC COEFFICIENTS 



Surge Hydrodynamic Coefficients : 



*;=o.o 


K= o-o 


X 

II 

o 

o 




X' w =-0.01743 


o 

o 

II 

v 8 

X 


X' u = -0.00282 


x 1=0.0 




*k=o.o 


X ' n =-0.00753 


x' wq =o.o 




<k=o.o 


X' ss =-0.01018 


X’s s =-0.01018 

°sp°sp 


^u=-° 01018 


X„=- 0.4024 







Sway Hydrodynamic Coefficients: 



II 

o 

o 


Y? =-0.03430 


o 

o 

II 

v* 


Yi =0.01241 


Y' =-0.00178 


II 

o 

o 


Y' =0.0 

wr 


Yi =0.01241 


II 

o 

o 


y; =o.oi i87 


Y' =-0.10700 


y;=o.o 


n 

o 

b 


Y' =0.0 

vw 







Heave Hydrodynamic Coefficients: 



Z' =-0.00253 


Z'. =-0.09340 


Z' w =-0.78440 


N 

•5 

II 

o 

o 


N 

II 

o 

o 


Z' =-0.07013 


N 

II 

o 

o 


Z'pr =0.0 


=-o.o2iio 


N 

II 

o 

o 


N 

II 

o 

o 


Z* p =-0.02 110 



Roll Hydrodynamic Coefficients: 



AT' =-0.00024 


<=0.0 


<p =0.0 


<=o.o 


K' p =-0.00540 


O O-O 


< 9 = 0.0 


<=0.0 


<=o.o 


K= 0-0 


< 9 =o.o 


K> O.o 
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Pitch Hydrodynamic Coefficients: 



M'. =-0.00625 


M'- =-0.00253 


m; =0.05 122 


a/; = 0.0 


M' q =-0.01530 


II 

o 

o 


“'3 >• 

II 

o 

o 


M' vp = 0.0 


M' &p =- 1.7664 


§ 

II 

o 

o 


II 

o 

o 


M; p = 1.3260 



Yaw Hydrodynamic Coefficients: 



at; =0.0 


=-0.00178 


a/; =o.o 


N' & =- 1.7663 


n; =-0.00047 


n; =o.o 


AC =0.0 


A^v =1-3259 


■3 > 

II 

o 

o 


N' r =-0.00390 


A' =-0.00769 


NL=0-0 


II 

O 

O 


<5 N 

II 

o 

o 







EQUATIONS OF MOTION 

The equations of motion for the NPS Phoenix AUV are presented in this section 
of the Appendix. These equations are based on the vector-matrix equation given in 
equation 2.87. The expanded equations use the assumption that Phoenix possesses both 
horizontal and vertical plane symmetry. 

Mass Matrix: 




m-X. 


0 . 


0 


0 


mz G 


-tny G 


0 


3 

1 


0 


-(mz c + y.) 


0 


rrvc G -Y , 


0 


0 


m ~ Z w 


my G 


~{ mx G + Zg) 


0 


0 


-(mZc+Ki) 


™}'g 






-(I« + Kr) 


mz G 


0 




~ I* 




-tyz 


-my G 


mx c -Y. 


0 






i 

l 
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Centrifugal/Coriolis Matrix: 



C(x) = 



0 


— mr 


mq 


m(y c q + z c r) 


- mx G q - Z^w 


mr 


0 


- mp 


-my c p + Z w w 


m(z c r + x c p) 


- mq 


mp 


0 


-mz G p ~y v v 


— mz c q + X ju 


-m(y c q+ z c r) 


mycP-Z^w 


mZcP + Y i v 


0 


- I yzq~ I xzP + ! z r ~ N t r 


mx c q + Z^w 


-m(z c r + x c p) 


mz c q~ X u u 


I y: q+ I Xi p - /./' + N-r 


0 


mx c r-y-v 


myc r + x v u 


-m(x c p + y c q) 




I zz r+I xyQ-txP+KpP 



Control Allocation Matrix: 












B(x)= 


UrX r6 br X S r S r ^br 


vp+ip +44 x 6 b p* bp 6 »p urX 


HH X *r*r S » 


U 4 X <*sp +tt H X *sp*sp d sp 


1 


i 


0 


U \ U \ Y *' 


0 


u\u\Y &r 


0 


0 


0 


0 


0 


U \ U \ Z »» 


0 


M Z * P 


0 


0 


i 


0 


0 


0 


0 


0 


0 


0 


0 


u\u\M Ap 


0 


“M m *p 


0 


0 


X tM 




0 




0 


-yu 


-y n 


0 



Control Input Vector: 



- mx c r + Y^v 
-my c r-X u u 
m(*cP + yc<l) 

-I~ r -l*y<l + lxP-KpP 
0 



0 0 0 

1 0 1 

0 1 0 

0 0 0 

0 -x„' 0 

X blt ® X sli 



u = [S br , 8 bp , S sr , 8 sp , F b , F rs , F bvl , F blt , F„ , F sll ] T 



Euler Angle Rates and Global Positions: 

X = U f +ucosy/cos6 + v[cos^sin 0sin 0 - sin y/cos<p] + wfaw y/sindcostp + sin y/sintp] 

Y = V f + u sin y/cosd + v[sin ^sin 0sin ^ + cos ^cos0] + vtfam y/sindcostp - cos y/sintp] 

Z = W f -usu\$ + vcos$sm<p + wcos6cos<p 

<p = p + qsin <p tan $ + rcos (ptan 6 

6 ' = qcostp- rsintp 

. _ (qsin<p + rcos <p) 
cose 
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Rotation Matrix: 



This rotation is based on equation 2.86. 








cos^cos# 

cos if/ sin # sin 0 - sin y/cos><p 
cos^sin 6cos(p + sin $/sin (p 



smy/cd 

sin ^sin #sin (p + cos y/cos(p 
sin ^sin#cos0-cos^sin0 



t 2 = 



1 

0 

0 



0 -sin # 

cos# sin # 

-sin0 cos#cos0 



-sin# 
cos# sin (p 
cos6cos(p 



Surge Motion Equation: 

m[u-vr + wq- x c (q 2 + r 2 ) + y G (pq -r) + z c (pr + q) ] - Z^wq + T-vr 

= X pp p 2 +X qq q 2 +X n r 2 +X pr pr+X a u + X wq wq + X vp vp + X vr vr 
+ uq(X qS J bp+ X q5 J sp ) + ur(X rS J br + X rSs S sr ) 

+ X m v 2 + X ww \v 2 + u\u\( X ss ( S 2 +S br )+ X 5sp5sp 8 2 sp + S bp ) 

-(W - B )sin6 + F ls + F rs - X res u\u\) 

.'lz£\t n u, +t a V, +T a W, +T n 0, +T„V, +T n W,) 
v S ) 
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Sway Motion Equation: 



m[ v + ur— wp + x G ( pq + r)-y G ( p 2 + r 2 ) + z G ( qr - p )] + Z-,wp + X -ur 
= Ypp + Y-r + Y pq pq + Y qr qr + Y,v + Y p up + Y r ur+ Y vq vq + Y wp wp + Y wr wr 



+ Y v uv + Y m vw + u\u\( Y &r S sr + Y &r S br ) 



P_ 

2 



UC dy h(x)\(v + xr^v + xr)] dx 

irtail 



+ (W -B )cos 0 sin (p + F blt + F sls 



r W -B '' 
v 8 J 




+ t 22 V, + t_ 3 W / +T,,U f +TvV, +T a w,} 



Heave Motion Equation: 

m[w-uq + vp + x G ( pr- q ) + y G ( qr + p )- z G ( p 2 + q 2 )] -Y-vp + X u ur 
= Z.q + Z pp p 2 + Z pr pr + Z n r 2 +Z^w + Z q uq + Z vp vp + Z vr vr 
+ Z vv v 2 + Z w uw + u\u\( Z &p S sp + Z^S^ ) 

- - J Jf C dz b( xjf w-xq )\( w - xq )] dx 

+ (W- B)cos<pcos0 + F bvl + F^ 

+ t»V, +i„w, +T n U, +T X V, +T»W,] 

\ 8 J 
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Roll Motion Equation: 



I x p+(I Z - I y )qr + I xy (pr - q)- 1 yz (q 2 - r 2 )- I Jpq+ r)- Z^vw + Y+vw 
- N-qr + M -qr + m[y G ( w-uq + vp)- z G (v + ur- wp)] 

= KpP + K-r + K pq pq+ K qr qr 

+ K+v + K p up + K r ur + K vq vq + K wp wp + K wr wr + K v uv + K m ,vw 
+ (y G W — y B B)cos<t>cosd - (z G W - z B B)cos0sin <p 



I y q + (I x ~ I z )pr ~ I xy (qr+ P)+ 1 yz (pq~ r) + 1 K (p 2 -r 2 ) 
-m[x G (w-uq+vp)-z G (u-vr+wq)] + Z w uw-X il uw + N i .pr-K p pr 
= M-q + M pp p 2 +M ^pr + M^r 2 + M ^w + M q uq + M Vp vp + M Vr vr 
+ M„v 2 +M w uw + u\u\(M &p S sp +M a , p S bp ) 



~(x G W -x B B)cos<pcosd-(z c W - z B B)smd-x bvl F bvt -x M F M 




Pitch Motion Equation: 




8 

r Wx c -Bx B 
K S 





J 
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Yaw Motion Equation: 



Ij + (Iy -I x)pq-lxy(p 2 - P 2 )- 1 yz(pr + q) + 1 Jqr - p) 

+ m[x G ( v + ur-wp)- y G ( u-vr + wq)] - F. wv + X-uv - M ij pq + K p pq 
= N p p + N,r + N pq pq + N qr qr + N,v + N p up + N r ur + N vq vq + N wp wp + N wr wr 
+ N v uv + N w vw + u\u\ N &r 8 sr + N &r S br ) 

- -j Jj C dy h( xj(v + xr\(v + xr)x\ dx 



+ (;c c lF — x B B)sin0cos0 + ( y G W - y B B)sind 

+ X bU F bU + X slt F sU ~ yis F ls ~ y rs F rs 



+ 



f Wy c -By, 
S 

v i 



fc.v, +t i 2 v, +i a w, +t h u, +r, 2 v, +T a w,\ 

\[t 2t U, +T n V f +T S W, +T,,U f +Tj t +T»W f ] 
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APPENDIX B. DOPPLER SENSORS 



This Appendix provides an overview of the two Doppler sensors used for control 
implementation in this dissertation; namely the SonTek® ADV Ocean and the RDI® 
Navigator DVL. 

The Sontek Advocean 

The SonTek ADVOcean (Acoustic Doppler Velocimeter Ocean Probe) is a 
versatile, high-precision instrument used to measure 3D water velocity in the most 
challenging applications, Figure B.l. The ADVOcean is designed for a wide range of 
environments including the surf zone, open ocean, rivers, lakes, and estuaries. 




Figure B . 1 ADV Ocean Probe 



The ADVOcean uses acoustic Doppler technology to measure 3D flow in a small 
sampling volume located a fixed distance (18-cm) from the probe. Figure B.2. The 
velocity range is programmable from ±5 to ±500 cm/s. Data can be acquired at sampling 
rates up to 25 Hz. 

With no zero offset, the ADVOcean can be used to measure flow velocities from 
less than 1 mm/s to over 5 m/s. The remote sampling volume, stability, and rapid 
response of the ADVOcean make it perfect for al 1 types of flow measurement: mean 
currents, waves, and turbulent flow parameters. 
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The ADVOcean consists of two elements: a probe and processor. The probe 
includes the acoustic sensors, receiver, and optional sensors in a submersible housing. It 
is connected to the processor using a custom shielded cable. 




Figure B.2 ADV Sampling Volume 

The ADVOcean processor operates from external DC power and outputs data 
using serial communication or a set of analog voltages. The processor can be operated 
from any PC-compatible computer or can be integrated with a variety of data acquisition 
systems. 

Standard Advocean 

The standard ADVOcean probe, Figure B.3, is designed for long-term 
deployments in hostile environments. The entire probe is constructed from 316 stainless 
steel, and protected from corrosion by a sacrificial zinc anode. With no moving parts, the 
ADVOcean has excellent resistance to biological fouling. For added protection, the entire 
probe, including the transducers, can be coated with anti-fouling paint. The probe is 
connected to the processor using an underwater mateable connector. 
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Figure B.3 Standard ADVOcean Probe 

For deep-water deployments, the ADVOcean can be rated for depths up to 2000 
m (the standard depth rating is 400 m). Deep-water ADVOcean systems are commonly 
used on autonomous underwater vehicles (AUVs) and remotely operated vehicles 
(ROVs) for detailed microstructure measurements. 

In any configuration, the ADVOcean probe is immune to zero drift and has no 
inherent minimum detectable velocity. The probe calibration can only change with 
physical damage to the system. No regular maintenance or re-calibration is needed. 
Advocean With Optional Sensors 

The ADVOcean probe can include a number of optional sensors to greatly expand 
its measurement capabilities. These include a compass and 2-axis tilt sensor allowing the 
ADVOcean to report velocity data in Earth (East-North-Up) coordinates; a pressure 
sensor for wave height estimation (PUVW) and surface-level measurements; and a 
temperature sensor for automatic sound speed compensation. 

ADVOcean probes with optional sensors use an expanded housing. Figure B.4, 
constructed from acetyl (Delrin), and have the same reliability and performance as the 
standard ADVOcean probe. 
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Figure B.4 ADVOcean Probe with Optional Sensors Housing 
Advocean Processor 

The ADVOcean processor can be housed in two different ways depending upon 
whether the processor will need to be submerged. The ADVOcean processor operates 
from DC power and is typically connected to a portable computer running SonTek's 
powerful data acquisition software. It can also be integrated with a variety of data 
acquisition systems using serial communication or the analog output voltages. 

Standard Features 



ADVOcean systems include the following standard features listed in Table B.l. 



ADVOcean Probe 


ADVOcean Processor 


• Factory Calibration (can only change 
with physical damage to the probe) 

• Programmable velocity range from 
±5 to ±500 cm/s 

• Submersible to 400 m 

• 10-m cable to processor (50-m 
max.) 


• Dual serial communication (RS-232 
standard, RS-422 for cable lengths to 
1500 m) 

• Four analog output voltages (3 
velocity, 1 signal strength) for 
integration with analog data 
acquisition systems 

• Hardware synchronization with 

external sensors using sync input and 
output 



Table B. 1 Standard ADVOcean Features 
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Options 

Several options are available for use with ADVOcean systems. The most 
common are the compass, pressure and temperature sensors. The compass and 2-axis tilt 
sensors allow the ADVOcean to report velocity data in Earth (East-North-Up) 
coordinates. The sensor has a built-in calibration feature to compensate for magnetic 
distortion. The user can easily re-orient the compass for up, down, or side-looking 
operation. Any ADVOcean with compass/tilt or pressure sensor includes a temperature 
sensor to compensate for changes in sound speed. Sound speed is used to convert 
Doppler-shift to water velocity. 

Advocean Performance Specifications 

The performance specifications of the ADVOcean are listed in Table B.2. More 
information may be found at the SonTek web site http://www.sontek.com. 



Performance Specifications 


General Operation 


Compass/Tilt Sensor 


• Acoustic frequency: 5 MHz 


• Resolution (heading, pitch, roll): 0.1° 


• Sampling rate: Programmable from 0.1 to 


• Accuracy (heading): ±2° 


25 Hz 


Accuracy (pitch, roll): ±1° 


• Sampling volume size: 2.0 cm3 


Temperature Sensor 


• Distance to sampling volume: 18 cm 


• Resolution: 0.01 °C 


• Minimum water depth: 20 cm 


• Accuracy: ±0.1 °C 


• Input power : 12-24 VDC 


Pressure Sensor 


Velocity Data 


• Available full-scale ranges: 10, 20, 50, 


• Range: Programmable to ±5, 20, 50, 200 


100, 200, and 500 m 


or 500 cm/s 


• Resolution: 0.00025 x (full scale) 


• Resolution: 0.01 cm/s 


• Accuracy: ±0.5 percent full scale 


• Accuracy: ±1 percent of measured 


Environmental 


velocity, ±0.25 cm/s 


• Operating temperature: -5° to 40°C 


• Random noise: Approximately 1 % of 


• Storage temperature: -10° to 50°C 


velocity range at 25 Hz; 


Dimensions 




• See Figure B.5 



Table B.2 ADVOcean Performance Specifications 
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Figure B. 5 ADV Ocean Dimensions 



The RDI Navigator DVL 

The WORKHORSE NAVIGATOR DVL is designed to provide rapid, precise 
velocity updates. Its small size, high accuracy, and low power consumption make it well 
suited to applications such as station keeping and sea bed surveys from underwater 
vehicles. The NAVIGATOR can be integrated with existing navigational systems 
(USBL, LBL and/or inertial). Its high-resolution velocity data provides a better focused 
picture of your vehicle's location and altitude. The NAVIGATOR is less than half the 
length and weight of standard broadband DVLs. Average power consumption of the 
WHN-1200 is only 8 watts. The NAVIGATOR is about 60 % of the price of standard 
broadband DVLs. The NAVIGATOR uses the patented Broadband technology. The 
WHN-1200 kHz has velocity accuracy of +/-0.2% +/- 0.2 cm/s. The NAVIGATOR 
measures velocity, altitude, heading, pitch/roll, and temperature. For a description of the 
operations of this sensor refer to [Gordon 1996] 
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WHN-1200 Specifications 



The specification for the 1.2 MHz Navigator are presented in Table B.3. Further 
specifications can be found at the RDI web site http://www.rdinstruments.com. 



Transducer and Pressure Case 


Actual Frequency 


1229 kHz 


Beamwidth 


1.2° 


Beam Angle (from vertical) 


30° 


Configuration 


4-beam-convex 


Housing & Transducer 


6061 aluminum 


Material 


External Connector 


7-pin low-profile 


Weight (in air) 


6.4 kg 


Weight (in water) 


4.2 kg 


Altitude 


Minimum 


0.3 m 


Maximum 


30 m 


Bottom Velocity 


Short Term Error 


0.3 cm/s 


(V = 1.0 m/s) 


Long term Error 


±0.2% ±0.2 cm/s 


Ping Rate 


1-10 Hz 



Table B.3 Navigator DVL Specifications 
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APPENDIX C. THE NFS PHOENIX AUV 



The physical layout of the NPS Phoenix AUV is shown in Figure Cl. Detailed 
description of the vehicle can be found in [Marco 1996] and [Brutzman 1997]. In 
addition an online description can be found on the NPS center for AUV Research web 
site at http://www.cs.nps navy.mil/research/auv/auvstats. 
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Drawn by D. Marco 1999 



Figure C. 1 Physical Layout of Phoenix AUV 

Prior to September 1996, Phoenix was used extensively as a test tank research 
vehicle. To transition the vehicle and the center to an ocean going operational unit 
required some significant upgrades in vehicle capabilities and center acquisitions. Table 
C.l summarizes the upgrades and acquisitions that were performed to allow Phoenix to 
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become deployable in the ocean. It also highlights the design features of the new vehicle, 
Figures C.la and C.lb, presently being outfitted at the center. 




Figure C.la Wire Frame Diagram of New NPS AUV 




Figure C.lb Solid Model of New NPS AUV 



188 





Pre-July 1996 


Present 


New Vehicle 


Propulsion Motors 


two 24v brushed DC 


two 1/4 hp, 24 v brushless 
DC (see figure C.2) 


two 48v DC, 1/4 hp, 
brushless 


Propellers 


3.5 in model submarine 
propellers (see figure C.3) 


7 in ducted propellers (in 
house designed) (see figure 
C.4) 


7 in ducted propellers (in 
house designed) 


Control Actuators 


eight control surfaces, (two 
fwd rudders, two aft rudders, 
and a pair of bow and stem 
planes), four 3 .5 in thrusters 
(two vertical fore and aft and 
two horizontal fore and aft) 
(see figure C.5) 


seven control surfaces, (one 
lower fwd rudders, two aft 
rudders, and a pair of bow 
and stem planes), four 3 .5 
in thrusters (see figure 
C.4a(two vertical fore and 
aft and two horizontal fore 
and aft) (see figure C.5) 


seven control surfaces, (one 
lower fwd rudders, two aft 
rudders, and a pair of bow 
and stem planes), four 6 in 
thrusters (two vertical fore 
and aft and two horizontal 
fore and aft) 


Vehicle Control 
Computer 


GESPAC Computer System 
running OS-9 real-time 
operating system 


GESPAC Computer System 
running OS-9 real-time 
operating system 


PC- 104 with a Pentium chip 
(166 MHz) running QNX 
(see figure C.6) 


Mission Control 
Computer 


Sun Solaris 


PC- 104 with a Pentium chip 
(90 MHz) running QNX (see 
figure C.6) 


PC- 104 with a Pentium chip 
(166 MHz) running QNX 
(see figure C.6) 


Electrical Power 
System 


two independent 24v lead 
acid battery systems 


two independent 24 v lead 
acid battery systems 


single 48v Absorbed Glass 
Mat (AGM) battery system 


Ballast System 


manual lead ballast 


manual lead ballast with 
syntactic foam for minor 


Variable ballast system 


Communication 

System 


Ethernet cable while in test 
tank 


Ethernet cable while in test 
tank, 900 MHz spread 
spectrum modem during 
missions (see figure C.7) 


Ethernet cable while in test 
tank, 900 MHz spread 
spectrum modem and u/w 
acoustics during missions 


Navigation System 


DR using water speed and 
vehicle heading 


EKF fusing Doppler and 
vehicle motion 


EKF fusing Doppler, GPS, 
LBL and vehicle motion 


Attitude Sensor 


three axis Mechanical Rate 
Gyro’s and a vertical gyro 


6DOF Solid State inertial 
sensing system (see figure 
C.8) 


6DOF Solid State inertial 
sensing system 



Table C.l Vehicle Upgrades and Acquisitions 
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Heading Reference 


Directional Gyro 


Directional Gyro with a 
PrecisionNav electronic 
compass backup (see figure 
C.9) 


Directional Gyro (possible a 
Honeywell ring laser gyro) 
with a PrecisionNav 
electronic compass backup 


Speed Reference 


turbo probe 


RDI DVL and SonTek ADV 
(see figures C. 1 0 and C. 1 1 ) 


RDI DVL and SonTek ADV 


Sensor Suite 


ST-725 Scanning Sonar, 
ST-1000 Imaging Sonar (see 
figure C.l 2) 


ST-725 Scanning Sonar, 
ST- 1000 Imaging Sonar, 
RDI DVL, SonTek ADV 


ST-725 Scanning Sonar, 
ST-1000 Imaging Sonar, 
RDI DVL, SonTek ADV, 
Altimeter , video camera 


Support equipment 


Transport cart (see figure 
C.l 3) 


Wells Cargo® Travel trailer 
outfitted with a mobile lab, 
12' inflatable boat and 
trailer, portable generator 
and computer workstation, 
vehicle shipping containers 
LXT acoustic tracking 
system (see figures C.l 4, 
C.l 5, C.16 and C.17) 


Wells Cargo® Travel trailer 
outfitted with a mobile lab, 
12' inflatable boat and 
trailer, portable generator 
and computer workstation, 
vehicle shipping containers 



Table C.l Vehicle Upgrades and Acquisitions (Cont.) 
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Figure C.2 Brushless DC Motors 





Figure C.3 Old 3-in. Props 
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Figure C.4 Present 7-in Ducted Props 




Figure C.4a Present 3.5-in Thrusters 
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Figure C.5 Phoenix Control Actuators 




Figure C.6 Mission Control Computer 
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Figure C.7 FreeWave Modem 
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figure C.9 PrecisionNav TCM2 Compass 




Figure C. 10 RDI Navigator DVL 
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Figure C. 1 1 ADVOcean Acoustic Doppler Velocimeter 




Figure C. 12 ST-725 and ST- 1000 Sonars 
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Figure C. 14 Mobile Lab Internals 
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Figure C. 15 Mobile Lab 




Figure C.16 Support Craft 
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Figure C.17 LXT Acoustic Tracking System 
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APPENDIX D. AUVFEST 98 



OVERVIEW 

The Naval Postgraduate School (NPS) participated in the NAVO sponsored 
AUVFEST for the first time, [Bunce 1998]. The lessons learned from this exercise were 
extremely valuable for future operational planning and vehicle development. First, the 
NPS Phoenix vehicle, a research platform, had never been deployed offshore from a 
research ship before. This task provided challenges that the four-man team failed to 
recognize beforehand, but was able to overcome. Like all other participants, the vehicle 
and support equipment had to be transported from their base of operations to Gulfport, 
MS. This movement of equipment was new to the Center for AUV Research 
(www.cs.nps.naw.mil/research/auv) , but enlightened the Center on the logistics involved 
with offshore operations. 

During the work-ups for this exercise, two significant hardware problems 
occurred. First, the Doppler unit originally integrated into Phoenix failed the week just 
prior to departure for Gulfport. This sensor was beyond repair, and a new Doppler 
needed to be purchased. A RDI Navigator DVL, a $25K unit, was purchased, through 
the Naval Supply System, and delivered to the Center in five days. The purchase of this 
unit in this time period was remarkable considering the government regulations that must 
be followed for a major purchase of this type. The vehicle, the support equipment, the 
new Doppler and all other sensors were shipped to Gulfport, installed, integrated into the 
vehicle control system and tested without faults in three days, a major accomplishment 
for a group that had never performed missions away from their base of operation. 

Secondly, during the weekend prior to shipboard load out, after the vehicle had 
been reassembled and all systems had been successfully tested and verified as 
performing, the SonTek Acoustic Doppler Velocimeter (ADV) was physically damaged. 
The damage to the ADV was not recognized until the after the vehicle had been loaded 
aboard the R/V GURE, and vessel was underway. The Center was able to contact the 
vendor and have a new unit shipped overnight to Gulfport where immediate configuration 
and installation was accomplished without error, affording Phoenix the ability to remain 
operational. 
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The objectives of the Center during AUVFEST ’98 were two-fold. First, NPS 
planned to use the sensors suite installed on Phoenix to demonstrate the ability of a 
moving platform to characterize the seaway. The theory and results from this section 
have been presented in Chapter VII and VIII, respectively. The second goal was to use 
the Phoenix’s forward-looking, sector-scanning sonar (Tritech ST725), to image water 
column targets and demonstrate the Phoenix’s onboard target identification algorithm. 
Navigation of Phoenix was accomplished using the RDI DVL together with a directional 
gyro heading reference in a deadreckoning filter. Details with regard to the above goals 
and vehicle missions are presented in the following sections. 

LOGISTICS 

The Phoenix AUV and all its associated support equipment was air freighted from 
Monterey, CA to Gulfport, MS via FedEx. This was a challenge for the Center since this 
was the first operational deployment. Custom shipping containers were purchased from 
Hardigg Industries, see Figure D.l, to hold the hull and nose fairing. These containers 
performed extremely well and protected the vehicle. The NPS shipping department 
provided packaging of the centers support equipment, see Figure D.2. Again due to the 
professional nature of the NPS employees, not a single item was damaged or turned up 
missing. 




Figure D.l Phoenix Shipping Containers 
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Figure D.2 Packaging of Support Equipment 



Upon arrival in Gulfport, the Center set-up shop in the NAVO BOATDET office, 
see Figure D.3. The vehicle was reassembled, including the newly purchased Doppler, 
bench tested and water tested, without any system degradations in two days. 




Figure D.3 Gulfport Temporary Work Space 



SHIPBOARD OPERATIONS 

The Phoenix, with its support equipment, was loaded aboard the Texas A&M 
University research ship, R/V GYRE , see Figure D.4. In addition to NPS, there were 
participants from Florida Atlantic University (FAU), Woods Hole Oceanographic 



203 



Institute (WHOI) and Massachusetts Institute of Technology (MIT) taking part in the 
demonstrations. 

i* 

* % 

* 1 * 

ii i 






Figure D.4 R/V Gyre at Anchor 

During the seven day exercise, the Phoenix conducted 27 separate open ocean 
runs, see Table D.l for a sample. This was a significant accomplishment seeing that the 
Phoenix was designed for tank testing some 10 years earlier. 

Operations from the ship were challenging. This was the first time Center 
personnel launched Phoenix from a crane, see Figure D.5. By the end of the second day 
of operations, the Phoenix crew was acting as if this was “old hat.” 




Figure D.5 Phoenix Being Launched from the Gyre 

The exercise was conducted in the Gulf of Mexico, south of Gulfport, MS. Figure 
D.6 displays the operating area. In this operating area, there were three different mission 
boxes. Each box had separate features so that the vehicles from the various participants 
could demonstrate their capabilities. 
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Figure D.6 AUVFEST Operations Area 



DIRECTIONAL WAVE INFORMATION 

Directional wave information was obtained during the AUVFEST missions. See 
Chapters VII and VIII for more information. 

AUTOMATIC SONAR TARGET RECOGNITION EXPERIMENTS 

One of the objectives of the AUVFEST experiments involved testing a recently 
developed sonar target recognition software module. The software is designed to process 
the sonar returns in real time and provide a reduced data set from the large amounts of 
data gathered. From the data, the centers of concentrations of high intensity sonar returns 
are identified and their locations stored during mission time, which requires no post 
processing. The location information can be then be used for post mission analysis or 
path re-planning to return to these areas for further study during the same mission. High 
intensity concentrations suggest underwater objects or structures, while areas of low 
intensity do not. Since the majority of the open ocean is clear, only a small amount of the 
data gathered is meaningful, and this approach greatly reduces the data storage 
requirements of the onboard computer system. The following will give a brief description 
of the sonar used, the identification algorithm, and finally the results from one of the 
missions. 
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Day/Data 
file name 


Run comments on Run 

Number 


Problems/ Fixes 


Sea State comments 


General info 












DAY 1 
1103.01 


1 Set up/balancing/vehicle 

was heavy by nose, -21 bs 
light. Time base run 4 mins 
(2 min out, 2 in back) at 3' 
depth. 


added buoyancy in nose / 
tow float/ reset fins/ 


light (1 foot seas, multi- 
directional, wind -10 mph) 


salinity 29.5 ppt at surface to 
34.0 ppt at bottom. Water 
temp -74 degrees. 


DAY 2 










1104.02 

1 104.03 

1 1 04.04 


1 TB run. At 090 true to go to 
GB. Initialize dg to ships 
head. Went to command 
but 180 off expected 
course. Mag compass 330 
degrees. 

2 repeat of first run. no 
changes to script. Same 
results as first run. 

3 Same run as first. 

Changed DG offset by 
2*pi. Results were no 
different than first run 


Fix sign error in DG offset 
calc. Due to sign error, we 
ended up with twice the 
ship's heading, which was 
100 degrees, error. 

Fix sign error in DG offset 
calc 

Fix sign error in DG offset 
calc 


sea state 1 -2 (2-3 feet, 4 foot 
swells). Seas from multi- 
directions. Swells from 115 
true. Wind driven from 025 
true. Wind 10-15 mph 


Taking ships heading/go to 
Gyre do a small WP run, a 
large WP run, then a run to 
area 2 buoy 


1104.05 


4 Same run as first. 

Attempted to give a 270 
heading command. Results 
were no different than first 
run 


Fix sign error in DG offset 
calc 






1 1 04.06 


5 Reinitialized vehicle 

headed at Gyre. New zero 
for DG. (-085 true) Time 
based run for 2 minutes 
with a commanded 
heading of 000. 


Behaved as commanded. 
Run was successful. 






1104.07 


6 Same run as previous. No 

reinitialization. 


Behaved as commanded. 
Run was successful. 






1104.08 


7 Time based run at -090 for 

five minutes. Heading 
away from Gyre. 


Behaved as commanded. 
Run was successful. 




Checked battery (24. 5v) and 
computer (22. 8v) voltages 
when mission was complete. 
Leak detectors at 1 .09. 


1104.09 


8 Time based run at 045 for 

seven minutes. Heading 
towards Gyre. 


Behaved as commanded. 
Run was successful. 




Comments from the chase 
boat was that we were doing 
-2.5 mph. 


1104.10 

1 


9 Time based run at 045 for 

two minutes. Heading 
towards the Gyre. 


Behaved as commanded. 
Run was successful. 






1 1 04.1 1 


1 0 small waypoint run for 90 
seconds heading away 
from the Gyre. Way points 
were (0,40,30), (-40,80,3), 
(-70,1 10,0) relative to 
initialized heading 


Behaved as commanded. 
Run was successful. 






1104.12 


1 1 long waypoint run. 

Attempted to start at (0,0) 
and returned to (0,0). 
Vehicle timed out trying to 
get to (100,-85) due to 
turning the long way. 


Run was moderately 
successful. 




Waypoints were (0,0), (50,0), 
(100,0), (150, -10), (190,-40), 
(200, -60), (180,-80), (140,- 
90), (100,-85), (50,-65), (20,- 
30), (0,0). With 6-meter 
diameter watch circle. 


1104.13 


12 long waypoint run. 

Attempted to start at (0,0) 
and returned to (0,0). 
Vehicle timed out trying to 
get to (100,-85) due to 
turning the long way. 
Attempted to increase 
watch circle. 


Run was moderately 
successful. Fix is to correct 
the heading command with 
an if statement to ensure 
that the vehicle always 
turns the shortest direction. 




Waypoints were (0,0), (50,0), 
(100,0), (150, -10), (190,-40), 
(200, -60), (180, -80), (140,- 
90), (100,-85), (50,-65), (20,- 
30), (0,0). With 10-meter 
diameter watch circle 



Table D.l Sample Phoenix Missions 
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SONAR HEAD GENERAL DESCRIPTION 

The NPS Phoenix is equipped with two mechanically scannable high frequency 
sonar heads built by Tritech Corp. One is a ST725 scanning sonar and the other is a 
ST 1000 profiler sonar. Each head can be scanned continuously through 360 degrees of 
rotation or swept through any defmed angular sector around the central axis of the unit. 
During normal operation, the head will ping, wait for return echoes to process, and then 

0 0 o 

rotate a specified angular width and repeat. Step widths of 0.9 , 1.8 , and 3.6 are 
computer selectable. 

All missions performed at AUVFEST ’98 used the ST725 which operates at a 

0 

frequency of 725 kilohertz and emits an acoustic beam 2.5 wide in the horizontal plane 

by 24° wide in the vertical plane. This device is described as a scanning sonar due to the 
nature of the range and intensity information returned for each ping. A scanning sonar 
operates by placing the intensities of the echoes from each ping into discrete segments of 
range called range bins. For this sonar, the number of range bins is nominally 128, but for 
operating ranges of 10 meters or less, the number of range bins is reduced to 64. The 
maximum operating range of the ST725 is 100 meters with a minimum operating range 
of 6 meters, and provides a resolution of (Maximum Range)/128 or (Maximum 
Range)/64, depending on the range setting used. The range setting used in the Gulf was 
20 meters, which gave a resolution of approximately 15 cm. 

In order to more clearly analyze the returns, the data can be thresholded to 
analyze only returns above a certain intensity level so that significant objects/structures 
can be seen, while other less significant entities (e. g. suspended particles in the water 
column, weak multi-path echoes, noise, etc.), are excluded. Combining thresholding with 
an appropriate power setting for the transducer, high quality results can be achieved. 
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SONAR CLUSTER IDENTIFICATION PROGRAM 



The identification algorithm is designed to recognize areas of contiguous high 
intensity sonar returns. Below is a section of a test case showing sonar scanlines that 
have been thresholded to record intensities above 10. 
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The algorithm records the centroid (X, Y pairs) of each cluster of high intensity 
returns, while ignoring noise or small concentrations such as the 14, 13 group shown 
above. Several parameters are selectable to tune the algorithm for target identification 
such as maximum cluster width, breadth, number of noncontiguous contacts, etc. The 
following presents the results of the identification algorithm from a run at AUVFEST. 
SONAR RESULTS FROM AUVFEST 

Since there were no submerged targets in the area where the Phoenix operated, the 
chase boat served as a suitable target. Figure D.7 shows the targets identified (chase 
boat) with the centroid of each marked with a cross-hair. The trajectory of the Phoenix is 
shown with a solid line while the sonar returns with an intensity above 10 are shown with 
asterisks. Figure D.8 shows the lower right target cluster with the centroid clearly 
identified by the algorithm. Since the AUVFEST results were very positive, further 
experiments will be conducted in Monterey Bay in the early summer of 1999. 




-150 -100 -50 0 50 



EAST (Meters) 

Figure D.7 Clusters Identified. 
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Figure D.8 Lower right cluster with centroid. 
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